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When using Wannier functions to study the electronic structure of multi-parameter Hamiltonians 

tf(fc.A) 

carrying a dependence on crystal momentum k and an additional periodic parameter A, 
one usually constructs several sets of Wannier functions for a set of values of A. We present the 
concept of higher dimensional Wannier functions (HDWFs), which provide a minimal and accurate 
description of the electronic structure of multi-parameter Hamiltonians based on a single set of HD¬ 
WFs. The obstacle of non-orthogonality of Bloch functions at different A is overcome by introducing 
an auxiliary real space, which is reciprocal to the parameter A. We derive a generalized interpo¬ 
lation scheme and emphasize the essential conceptual and computational simplifications in using 
the formalism, for instance, in the evaluation of linear response coefficients. We further implement 
the necessary machinery to construct HDWFs from ab initio within the full-potential linearized 
augmented plane-wave method (FLAPW). We apply our implementation to accurately interpolate 
the Hamiltonian of a one-dimensional magnetic chain of Mn atoms in two important cases of A: 
(i) the spin-spiral vector q, and (ii) the direction of the ferromagnetic magnetization m. Using the 
generalized interpolation of the energy, we extract the corresponding values of magneto-crystalline 
anisotropy energy, Heisenberg exchange constants, and spin stiffness, which compare very well with 
the values obtained from direct first principles calculations. For toy models we demonstrate that 
the method of HDWFs can also be used in applications such as the virtual crystal approximation, 
ferroelectric polarization and spin torques. 

PACS numbers: 71.15.-m, 75.75.-c, 77.84.-s 


I. INTRODUCTION 

Maximally localized Wannier functions (MLWFs) have 
become a widely applied tool in electronic structure cal¬ 
culations p]. Defined as discrete Fourier transformations 
of Bloch states 'll km with respect to crystal momentum 
k, the MLWFs 


W Rn (r) = 2- £ e~ ik R U^km(r) (1) 

km 

are labeled by the direct lattice vector R and the orbital 
index n. Unitary gauge transformations U^ as well as 
the number of fc-points, IV*., enter Eq. Q. These orbitals 
allow for an efficient but remarkably accurate Wannier 
interpolation of any single-particle operator such as the 
Hamiltonian . The Wannier interpolation is in par¬ 
ticular fruitful in the calculation of linear response co¬ 
efficients such as the anomalous Hall conductivity, and 
various Fermi surface properties, which require a fine k- 
mesh for Brillouin zone (BZ) integration [2H3]. 

It is sometimes necessary to consider a family FU fc,A ) of 
Hamiltonians, where A is an additional parameter. In the 
problem of ferroelectric polarization, for instance, the pa¬ 
rameter A indicates relative displacements of the crystal 
sublattices mu- In magnetic systems with non-collinear 
spin-spiral texture, the additional parameter A can be 
identified with the spin-spiral vector q. The related en¬ 
ergy E(q ) serves to determine Heisenberg exchange con¬ 
stants [8]. Frequently, the ferromagnetic magnetization 
direction rh plays the role of A. Such a situation is met 


in the study of the magneto-crystalline anisotropy energy 
(MAE), which is the magnitude of the variation of the 
energy E(rh) with fh. The dependencies of crystal vol¬ 
ume, current-induced torques [SHE], and the conductiv¬ 
ity tensor on the magnetization direction provide similar 
examples. 

Notably, the anomalous Hall effect (AHE) can exhibit 
an anisotropy with respect to rh [13] . For a given magne¬ 
tization direction, the corresponding value of the anoma¬ 
lous Hall conductivity is obtained from Wannier inter¬ 
polation in crystal momentum k. The evaluation on a 
dense rh-mesh requires accordingly the construction of a 
huge amount of MLWFs - one set of MLWFs for each 
rh. Consequently, the accurate calculation of the AHE 
anisotropy is a rather time-consuming task. 

The computation of such linear response quantities 
would benefit in particular from an interpolation tech¬ 
nique based on functions which provide efficient access 
to the multi-parameter Hamiltonian 7U fcA ^ of the sys¬ 
tem. For this purpose, the definition of MLWFs, Eq. 0 , 
has to be generalized. We introduce higher dimensional 
Wannier functions (HDWFs) as Fourier transformations 
of states with respect to both crystal momentum 

k and the additional parameter A: 




1 1 


e- ik R e~ lX S 

k\m 


x^„ A) $ fcAm (r,£). 


( 2 ) 


Here, ^ fc,A ) denotes a unitary matrix and the new ad¬ 
ditional index E of the HDWFs is conjugate to A like 
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R is conjugate to fc. Further, N\ is the number of 
A-points and £ refers to an auxiliary space variable. 
In the presence of an additional parameter, the usual 
Bloch states ^kXm are typically not orthogonal, i.e., 
(^feAnl^fc'A'm) 9^ S k k'5\\'5 nm . To overcome this obsta¬ 
cle and establish the transformation Eq. ([2]), the intro¬ 
duction of an auxiliary space £ is crucial. The auxiliary 
space is reciprocal to A like real space is reciprocal to 
crystal momentum k. Then, the role of Bloch states is 
taken in Eq. ([2]) by orthogonal states $ k xm in the com¬ 
bined space of r and £. 

Within an energy window of interest, the family H^ k ' X ^ 
of Hamiltonians can be interpolated using HDWFs. After 
constructing HDWFs from a coarse (fc, A)-mesh, we store 
the information on the multi-parameter Hamiltonian in 
hopping elements H nm (R. 3 ) of HDWFs. Finally, we 
can obtain H^ k ' X ^ on a much denser (fc, A)-mesh by an 
inverse Fourier transformation of these hoppings. 

Several applications where such an approach would be 
very fruitful come to mind, of which we mention explic¬ 
itly the following ones: (i) The evaluation of the AHE 
anisotropy would be simplified by performing a general¬ 
ized Wannier interpolation for the situation with A = rh; 

(ii) Heisenberg exchange constants would be accessible 
using generalizations of MLWFs in an interpolation of 
E(q) with respect to the spin-spiral parameter A = q: 

(iii) Mixed Berry curvatures in real and momentum space 
have been found to be quantitatively important in mate¬ 
rials like MnSi, where they support the formation of non¬ 
trivial magnetic textures M- An accurate interpolation 
of the multi-parameter Hamiltonian could be employed 
in the study of the contributions of mixed Berry cur¬ 
vatures to the Hall effects. Thus, HDWFs would prove 
useful in the topological characterization of complex mag¬ 
netic structures; (iv) Such functions could provide an al¬ 
ternative means of calculating forces in first principles 
methods where atomic displacements are described by 
A; (v) Eventually, the framework could allow the treat¬ 
ment of alloys like Fe^Coi^ or even Bi^Sbi-a, within the 
virtual crystal approximation (VCA), with concentration 
x as parameter A. 

In this work, we present the formalism of higher dimen¬ 
sional Wannier functions (HDWFs) given by Eq. <§• The 
problem of non-orthogonality of Bloch states is solved by 
the introduction of an auxiliary space reciprocal to the 
A-space. Based on HDWFs, we establish a generalized 
interpolation scheme which provides efficient but accu¬ 
rate access to the multi-parameter Hamiltonian H^ k,x ^ 
for any desired value of (fc, A). The necessary machin¬ 
ery for an ab initio construction of HDWFs is imple¬ 
mented within the FLAPW method to treat consistently 
multi-parameter Hamiltonians of realistic systems. As 
proof of principle, we consider the electronic structure 
of a linear equidistant chain of Mn atoms as a function 
of (i) the spin-spiral vector q , and (ii) as a function of 
the ferromagnetic magnetization direction rh. Using the 


method of HDWFs, we achieve the generalized interpo¬ 
lation of the first principles Hamiltonian family ijW'd 
and H^ k ’ m \ which allows for a precise determination 
of Heisenberg exchange parameters, spin stiffness and 
magneto-crystalline anisotropy energy. Within toy mod¬ 
els we investigate further promising applications of the 
formalism such as VCA, current-induced torques, and 
ferroelectric polarization. 

The paper is structured as follows. We begin with 
a concise review of MLWFs and the Wannier interpo¬ 
lation in Sec. [TTJ In Sec. m we introduce the formal¬ 
ism of HDWFs and set up the interpolation technique of 
multi-parameter Hamiltonians. We describe the imple¬ 
mentation for constructing HDWFs from ab initio within 
the FLAPW method in Sec. |IV[ and present the applica¬ 
tion of HDWFs to calculating Heisenberg exchange con¬ 
stants and MAE of a Mn chain in Sec. |IV| and|Vl respec¬ 
tively. In Sec. VI we discuss applications of HDWFs 
for VCA, ferroelectric polarization, and current-induced 
torques based on toy models. Finally, we conclude this 
work with a summary. 


II. REVIEW OF MLWFs 

In contrast to the oscillatory and delocalized Bloch 
states ’Ffc n , Wannier functions (WFs) provide a more 
intuitive insight into the nature of crystal bonding D3K 
ITS] and the underlying physical processes due to their 
real-space localization. The benefit of reformulating the 
electronic structure problem in terms of WFs is widely 
exploited in formal developments such as effective model 
Hamiltonian construction for the study of strongly corre¬ 
lated systems dum. Further, the centers of WFs play 
a fundamental role in the modern theory of ferroelectric 
polarization m- 

In the definition of MLWFs, Eq. 0. the unitary ma¬ 
trices U are chosen to maximize the real-space local¬ 
ization whereby the resulting orbitals are uniquely deter¬ 
mined. One approach to obtain such a gauge thus lies in 
minimizing the spatial extent Q of the WFs: 

n = £ ((Wo n |r 2 |Wo„> - (W 0n \r\W 0n } 2 ) . (3) 

n 

An algorithm for the spread minimization of WFs was 
proposed first for the case of isolated groups of energy 
bands [22] but soon generalized to treat entangled bands 
as well [23]. The corresponding WANNIER90 implementa¬ 
tion requires as an input two quantities [241 . First, the 
overlaps 


M mn^ = ( u k m\uk+bn ) (4) 

of the periodic parts Uk m = ^~ lk r ^km of the Bloch 
states at neighboring crystal momenta fc and fc + b have 
to be provided since they determine centers and spreads 
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of MLWFs. Second, the projections Ami = km\gn) of 
the Bloch functions onto localized trial orbitals g n serve 
as a starting point for the iterative minimization process, 
which results at the end in MLWFs. 

The Wannier interpolation is performed within a cer¬ 
tain energy window spanned by the MLWFs [3]. For this 
purpose, matrix elements of the single-particle Hamilto¬ 
nian H between such functions have to be calculated: 


H nm {R) = (W 0 n\H\W Rm ) 


1 


J2e~ lkR 



(5) 


where £ kn > stand for the ab initio band energies computed 
on a coarse fc-mesh of N k points. Importantly, MLWFs 
are orthonormal such that {Wu n \Wn’ m ) = SRR'S nm , 
which follows from the orthogonality of the Bloch states 
('Ffcnl'hfc'm) = N k 5 kk '5 nm and Eq. 0. Because of the 
localization of MLWFs, the matrix elements H nm (R) de¬ 
cay rapidly with increasing distance |i?|. The electronic 
band structure can be accessed accurately on a much 
finer interpolation mesh of fc-points using the hopping 
elements, Eq. 0. By an inverse Fourier transformation 
the interpolated Hamiltonian H is obtained for every 
desired fc-point, even if this point is not contained in the 
coarse mesh of N k points used for constructing MLWFs: 


Yl e ik RH nm(R) ■ ( 6 ) 

R 


Here, as marked with a dash, the summation is truncated 
keeping in mind the rapid decay of the hopping elements 
H nm (R). Eventually, the interpolated Hamiltonian is 
diagonalized using unitary matrices V^: 


y( fc ) 


— £ k .r)b 


kn^nm • 


(7) 


Thus, the Wannier interpolation grants efficient access to 
the band structure £ kn for any k. Key properties neces¬ 
sary for this interpolation scheme to work are orthonor¬ 
mality as well as real-space localization of the MLWFs. 


III. EXTENSION OF THE FORMALISM 
A. Orthogonality problem 


where £ k \ n are the band energies. Since the Hamiltoni¬ 
ans H A) and are generally independent, the eigen¬ 

states at different values of A are not necessarily orthog¬ 
onal, i.e., 


| A'm) 9^- ^kk' ^XX' $nm ■ (9) 

Only at fixed parameter A the orthogonality with re¬ 
spect to crystal momentum is always present such that 
k\n\^k'Xm) = N k 5 kk '5 nm ■ As a consequence, discrete 
Fourier transformations of these Bloch states with re¬ 
spect to k and A do not lead to orthonormal WFs. On 
the one hand, nonorthogonal WFs can be defined |J5] and 
can even be advantageous due to a stronger real-space lo¬ 
calization [2S]. On the other hand, in our case already the 
eigenstates are nonorthogonal for A /A', leading to ad¬ 
ditional complications. In particular, when trying to gen¬ 
eralize Eq. 0 for the case of these nonorthogonal WFs, 
we formally encounter matrix elements fitk\n\H\^ k \'m) 
the handling of which is not obvious for A ^ A'. 


B. Solution to the orthogonality problem 

1. Introduction of an auxiliary space 

To obtain well-localized orthonormal HDWFs, we in¬ 
troduce an auxiliary space £ as the reciprocal of the A- 
space. Instead of taking the Bloch states 'S’ k \ n (T') in the 
construction of HDWFs, we consider orthogonal states 
$k\n(r,$ l ) in the composite space (r,£). We define such 
states as the products of the physical Bloch states and 
auxiliary orbitals Ca(£) : 

®k\ n {r,£) = *kXn(r)UZ) ■ ( 10 ) 

The crucial orthogonality of the product states <FfcAn is 
enforced by choosing (Ca|Ca') = N x 6\\>. 


2. Choice of the auxiliary orbital 

When constructing the auxiliary orbital, we consider a 
translationally invariant potential in the auxiliary space 
as schematically shown in Fig. [T] The auxiliary orbital 
Ca(£) is chosen to be the lowest energy eigenstate of an 
according lattice periodic Hamiltonian H: 


In the presence of an additional periodic variable A, 
the system under consideration is described by a family 
of Hamiltonians, where each member HA) represents the 
system at a given value of A. If we assume that is 
lattice periodic at each A, the eigenstates of H ix > are 
Bloch states 'h k \ n carrying a dependence on A: 

H (X) {r)'H k \ n (r) = £k\n^k\n(r) , (8) 


m)C a(0 = £aCa(£), (11) 

where £\ represents the lowest energy band in A-space 
associated with the Hamiltonian H. The regular lattice 
in the auxiliary space is modeled using a series of po¬ 
tential wells of depth V 0 as depicted in Fig. [l] Because 
the extension to higher dimensions is straightforward, we 
only discuss the one-dimensional case described by the 
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single-particle Hamiltonian 


3 

where m is the electron mass and h is Planck’s constant. 
To simplify notation, we introduced the function 

Ol (0 = 0(£ — 3j+ 6/2) -G(€-3j- 6/2) , (13) 

which cuts out the well region of width b centered around 
the position 3j with the Heaviside step function 0(£). 
Here, the coordinate 2 j = ja is defined by the lattice 
constant a measured along the £ axis and an integer j. 

The most convenient and natural choice of the auxil¬ 
iary orbital is that of a Bloch wave: 

CaCO = e lA 'WD , (14) 


where px{£) is a ^-periodic function normalized to the 
unit cell in the auxiliary space: (px\px') = ^xx'- It 
follows that the auxiliary orbitals are orthogonal, i.e., 
(Ca|Ca') = N\5 \\>, where the integration is performed in 
a supercell of N\ unit cells in the auxiliary space. 

We can solve the Schrodinger equation to the one¬ 
dimensional Hamiltonian Eq. (12) for £a numerically us¬ 
ing a plane-wave basis, with the potential depth Vq cho¬ 
sen to strongly suppress the tunneling between different 
wells. Alternatively, we can also arrive analytically at the 
expression for the lowest energy eigenstate to Eq. (12) in 
the deep-well limit Vq —i oo by starting from a single-well 
solution: 


w(0 = 


0 


2cosy , if |^| <6/2 
, else 


(15) 


The auxiliary orbital is found as an inverse Fourier trans¬ 
formation of the Wannier-like function w with respect to 
the positions 2 ? . In accordance with the Bloch theorem, 


Eq. (14), the lattice periodic part p\ assumes the form 



Figure 1. Scalar potential landscape (red solid line) of the one¬ 
dimensional lattice defined as series of finite potential wells of 
depth Vo. The well width is b, and a stands for the lattice 
constant. A dashed line indicates the £ axis. 


3. Product states and composite Hamiltonian 


Essentially, the definition of HDWFs, Eq. ([2]), is based 
on the products §kXn('r,£) of Bloch states and the aux¬ 
iliary orbitals discussed above. Exploiting Eq. (fl4|), we 


can rewrite the product states of Eq. (101 as 


$k\n{r, £) = e lk r e lX ^ip k \n(r, £), 


where 


<Pk\n(r, £) = u k \ n (r)px($,) 


(18) 


(19) 


are lattice periodic. Such product states are orthogonal 
also in A, i.e., kXn\*& k'x'm} — -^k-^x^kk'^xx'^nrm and 
they are periodic with respect to both k and A. 

The question arises to which Hamiltonian PL the prod¬ 


uct states 'FfeA n(r,|), Eq. (18), are eigenstates in the 
composite space (r,£). Since the eigenstates have the 
product shape, the sought Hamiltonian decomposes into 
two additive contributions. If we denote by H the single¬ 
particle Hamiltonian to which the auxiliary orbital (a 
is an eigenstate (see Eq. ([TT])), the Hamiltonian of the 
composite system is given by 




( 20 ) 


Pv(C)=e w(£ - 2 j). 

j 


(16) 


Here, the Hamiltonian H , which is independent of the 
parameter A, can be written in the form 


Eventually, overlaps between periodic parts p\ to differ¬ 
ent parameter values A and A + t are important ingredi¬ 
ents for constructing HDWFs. Such overlaps read 


T 4 ), 

(17) 

where r plays a similar role like b in Eq. Q and the 
£ integration is performed in one unit cell in auxiliary 
space. From the above Taylor expansion it follows that 
{px\d\P\) = 0 and {px\d 2 x p\) = (6 - 7t 2 )6 2 /(127t 2 ). 


{p\\p\+r) = 


87t 2 sin (rb/2) 
47 T 2 rb — r 3 b 3 


= i+ S <6 -’ r2)+o( 


H{r) = J H (x \r)6(X- A) dA, (21) 

where A'FfcAn = ATfcAn- When acting with this Hamilto¬ 
nian on a specific Bloch state 'i’kXn, the delta function se¬ 
lects the Hamiltonian which corresponds to the spe¬ 
cific parameter value of the Bloch state. It thus follows 
that H'i’kXn = ZkXn'&kXn in line with Eq. |8|. There¬ 
fore, the product states satisfy the Schrodinger equation 

'H(r,£)$kXn(r,£) = (EkXn + £x) $kXn{r ,£), (22) 






























5 


where £\ represents the energy band in A-space associ¬ 
ated with the auxiliary orbital Ca- 


According to Eq. (22), the eigenvalues of the compos¬ 


ite Hamiltonian T-L differ from the ab initio band energies, 
which we would like to interpolate. To achieve the iden¬ 
tity between the two sets of eigenvalues, we study the 
deep-well limit for the Hamiltonian H. In this case, the 
energy level £\ becomes independent of A. As a conse¬ 


quence, £\ in Eq. (22) can be set to zero, so that 


H(r, £)&k\n(r, £) = £k\n$k\n(r, £). (23) 


Therefore, the generalized interpolation in k and A of the 
band structure of the composite Hamiltonian T~L grants 
access to the interpolated band structure of the physical 
Hamiltonian H of interest. 


C. Higher dimensional Wannier functions 
(HDWFs) 

1. Definition 


Discrete Fourier transformations of the product states 
with respect to k and A define HDWFs in close 
analogy to MLWFs. We repeat here Eq. <© as one of the 
main results of this work: 


W RSn {r,£) 


1 1 
Nk Nx 


Y e~ ik R e~ iX S 

kXm 


xU£ n x ^kXm(r,Z). 


(24) 


HDWFs are labeled by an orbital index n, the direct lat¬ 
tice vector i?, and an additional lattice vector S. S is 
conjugate to A like the direct lattice vector R is conjugate 
to the crystal momentum k. Unitary gauge transforma¬ 
tions U( k ' X ^ control the localization of HDWFs, and Nk 
and N\ stand for the number of grid points in fc-space 
and A-space, respectively. Due to the orthogonality of 
the product states <&kXmi the orbitals WRs n (r,£) are 
orthonormal, i.e., (WRs. n \W R "=' m ) = $rr' fes' S nm . 


2. Centers and spreads 


A first physical interpretation of the functions Wrsu 
defined by Eq. (24) is provided by the expressions for the 
centers of HDWFs in r and £. The centers of HDWFs in 
real space r can be directly related to the BZ sum of the 
Berry connection in crystal momentum space: 


which is easily derived using the Bloch-like periodic parts 

iPkXn = 7L Idm.n •pkXtn and UkXn = TC . 


, ,(fc. A) 

, Otmn UkXmi 


respectively. Equation (25) is the generalization of the 
expression for centers of MLWFs. To obtain the ^-centers 
of HDWFs, we start from the definition Eq. (24) and 
write down the expectation value of the auxiliary position 
operator in the basis of HDWFs: 


(W00n|€|W00n) = ^(yfcAnlVAly, 


k\n/ 


k X 


N k N x 


X ({UkXn\Vx\u k Xn) 


fcA 


(pa|V a |pa)) • 


However, the second term (pa|Va|pa) vanishes in the 
deep-well limit (see e.g. Eq. ©). Accordingly, the cen¬ 
ters of HDWFs in the auxiliary space £ are given by the 
BZ sum of the Berry connections in A-space: 

(WoOnlSIIUoOn) = AT AT X (UkXn | V A \u k Xn) i (27) 


k\ 


which are independent of the auxiliary orbitals Ca but 
determined solely by the Bloch-like periodic parts. 

Likewise, the expectation values for the squared posi¬ 
tion operators r 2 and £ 2 evaluate to 


{WoQn\r 2 \Woon) — 


-1 

-1 


NkNx fr 

k X 


kXn\ VfclpfcAra) 
fcA 


(28) 


(f^fcAn l^fc l^fcAn) ? 


and 


{W / 00 n|£ 2 | Woon) — 


-1 

XX 

-1 

NkNx 


{fikXn | ^ xl^fcAn) 


fcA 


^ ^ ((fifcAnl^ 


A l^fcAn) 


(29) 


fcA 


+ (pa|V||pa». 


While the ^-center is independent of the auxiliary orbital, 
the expectation value of £ 2 contains explicitly a con¬ 
tribution from the integral (pa|V||pa). Together with 
Eq. (25) and Eq. (27) for the centers, the above expres¬ 
sions can be used to calculate the spread £1 of HDWFs 
in the combined space of r and £: 


& — ((WoOn|7 2 |IUoOn) ~ (W'ho„ |f |Uoon)^ 


(30) 


To simplify notation, we introduced the generalized po¬ 
sition operator f = (r,£). 


(fUoOn | ^ | f f^OOn) 


AT jy ^ , {fPkXn l^fc I AfcAn) 

k X kX 

{UkXn | ^fc |fifcAn) j 


( 25 ) 


3. Maximal localization 

As discussed in Sec. [TIJ the constraint of minimal 
spread Cl uniquely defines the MLWFs up to a global 
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phase factor. Similarly, the unitary matrix IA in the 
definition of HDWFs, Eq. (24), is determined from the 
condition that the orbitals W R3n (r,£) should exhibit a 
minimal spread Q in the space of r and The resulting 
HDWFs are unique up to a global phase factor. 

Usually, the maximal localization procedure is per¬ 
formed in the three-dimensional real space r. For con¬ 
structing HDWFs, we need to consider in addition the 
auxiliary space £. To minimize the spread f2, Eq. (30), 
in the composite space of r and £, we thus extend the 
WANNIER90 program. Then, centers of HDWFs possess 
additional coordinates, Eq. (27), owing to the auxiliary 
space. A higher dimensional but block diagonal Bravais 
matrix is employed to define the composite direct lattice 
in (r, £)-space: 


A = 


A i 0 

o a 2 


(31) 


where A\ is the usual 3x3 Bravais matrix of the crys¬ 
tal and the rank of A 2 is given by the dimension of £. 
Associated with the direct lattice is a composite recipro¬ 
cal lattice in (fe, A)-space combining crystal momentum 
and the additional parameter. Both fc and A are cho¬ 
sen to form individual Monkhorst-Pack grids. According 


to the equations in Sec. Ill C 2 we can exploit the same 
finite-difference expressions as in Ref. [22] to evaluate the 
spread H. However, the role of the usual periodic parts 
u kn( r ) i s now taken by their higher dimensional analogs 
Vk\n{f ,£)> Eq. (19). We need to set up all necessary 
neighbors (k + b*,, A + b\) of a point (fc, A) to apply the 
finite-difference formulas. Here, bk and b\ connect the 
two reciprocal points. In general, we choose the Bravais 
matrix such that only those neighbors need to be consid¬ 
ered where either bk = 0 or b x = 0. 

The overlaps (pkXm\pk+b k x+b x n) of the periodic parts 
at neighboring points in the (fc, A)-space serve to calcu¬ 
late centers and spreads of HDWFs. As we choose the 
directions of k and A to be orthogonal in the composite 
reciprocal lattice, the overlap matrix consists of the two 
contributions 


A) = (VkXmWk+bXn) (32) 

M£V{k) = {VkXmWkX+bn) (33) 


which is discussed in Ref. m- The evaluation of the 
overlaps in Eq. (351 requires the integrals (px\px+b) in 
addition to the overlaps Mmn\k) = ( u k x m \ukx+bn ) 
between the periodic parts of Bloch states. Details on 
the implementation of M.mn\k) within the FLAPW 
method are given in the Appendices [A]fD[ for various re¬ 
alizations of A. 

In addition, the projections of Bloch states onto local¬ 
ized trial orbitals g n (r) are replaced by the projections 


= (<S>kXm\Pn) (36) 

of the product states onto functions Pn(»b£) localized 
in (r, £)-space. Such projections are the starting point 
for the spread minimization of HDWFs. Exploiting the 
product shape of the states Eq. ( [14] ), we obtain 

A& x) = (*kx m \g n )(tx\h), (37) 

with the ansatz p r ,(r, £) = < 7 n(r)/i(£). Thus, the projec¬ 
tions onto localized trial functions factorize into the usual 
projections of Bloch states and the auxiliary projection 
(Cx\h). In Appendix [Aj the construction of the usual 
projections (’ %kXm\9n ) within FLAPW is discussed. 


D. Generalized Wannier interpolation 


Within the energy window spanned by HDWFs, the 
multi-parameter Hamiltonian H^ k ’ X ^ can be interpolated 
in fe and A. As a starting point for the interpolation 
scheme, the matrix elements of the Hamiltonian T~L in 
the basis of HDWFs have to be calculated: 


H nm (R, S) = (W 0 on\n\W RSm ) 


1 

N k N x 


Y e~ ik R e~ iX ' s 

k\n' 


X 



£k\n' 


j /(k A) 
U n'm ' 


(38) 


These generalized hoppings are rapidly decaying with dis¬ 
tance and, further, depend only on the distance vectors 
R! — R and E' — E: 


depending on whether overlaps at neighboring fc-points 
or neighboring A-points are concerned. The product 
shape of the periodic parts ipkXm, Eq. (19), allows further 
simplifications: 


= (UkXm\Uk+b Xn) , (34) 

Mmn ^ (k) = (UkXm\ u k X+bn) (Px\PX+b) . . 

= M^ ] (k){px\px+b) ■ 


(w RSn \n\w RI3 r m ) = (w 0 on\n\w R '_ R3 '_ 3m ). ( 39 ) 

The hoppings H nm (R , E) converge quickly with the num¬ 
ber Nk x Nx of mesh points. They can therefore be 
constructed using a coarse Nk x N\ mesh. By an in¬ 
verse Fourier transformation one obtains the interpolated 
U( fc,A ) for every desired point (fe, A), even if this point is 
not contained in the coarse Nk x Nx mesh used for the 
construction of the HDWFs: 


The implementation of Eq. (34) within the FLAPW 
method is analogous to that of the usual overlaps, Eq. Q, 


H^ m X) = E' e lk R e lX S H nm (R , E) . (40) 

RS 
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As the HDWFs are strongly localized after minimizing 
their spread, the summation can be truncated. Indicated 
by the dashed symbol, only non-negligible hoppings are 
taken into account. Finally, the interpolated bands £ k \ n 
are obtained by diagonalizing the interpolated Hamilto¬ 
nian 


^ jj(k, A) y(k,\) 


kXnv 


(41) 


IV. APPLICATION TO SPIN SPIRALS IN A 
CHAIN OF Mn ATOMS 


where + kqn (r) and u k (r) are lattice periodic functions. 
Using the latter ansatz allows us to avoid computation¬ 
ally demanding first principles calculations of large su¬ 
percells and to perform all the calculations in a unit cell 
of one Mn atom. The Bloch states, Eq. (45), can be 
chosen to obey the periodic gauge in k and q simul¬ 
taneously. However, we emphasize that due to the in¬ 
dependent phases in Eq. (45), which arise from spin-1/2 
rotation matrices, the period associated with the q-mesh 
is enhanced by an overall factor of two (recall that spin- 
1/2 acquires a Berry phase of 7r upon rotating by 360°). 
Consequently, for the construction of HDWFs we have to 
uniformly sample the range [0,47r/a) of q-values. 


A. Heisenberg model and generalized Bloch 
theorem 


B. Symmetry 


As an application of the generalized Wannier interpo¬ 
lation to realistic systems, we study a one-dimensional 
magnetic chain of Mn atoms oriented along the z di¬ 
rection and extract Heisenberg exchange constants from 
HDWFs. The Heisenberg model is defined as 

H = — JijSi ■ Sj , (42) 

ij 


where the Heisenberg exchange constants Jij mediate the 
exchange interaction between the normalized moments 
S t and Sj located at the sites i and j, respectively. In 
case of the magnetic monatomic chain, the most general 
solution to Eq. (42) is the non-collinear (flat) spin-spiral 
state 


S n = (cos naq, sin naq, 0), (43) 


which is characterized by the spin-spiral wave vector q = 
qe z . Here, a is the lattice constant along the chain axis. 
If we exploit translational invariance Jij = Jo|j—»| and 
Eq. (43), the energy of the system assumes the form 


E(q) = -2 ^2 J 0n cos (naq). 


(44) 


Expanding the energy in the vicinity of the ferromagnetic 
state (q = 0), we can define the spin stiffness A of the 
magnetic chain through E(q —> 0) « E(Q) +Aq 2 . In order 
to access efficiently the Heisenberg exchange constants 
Jo n as well as the spin stiffness A, we treat the spin-spiral 
vector q as an additional variable of a multi-parameter 
Hamiltonian H^ k q \ 

Without spin-orbit interaction we can make use of the 
so-called generalized Bloch theorem, which dictates a 
specific Bloch-like shape of the spinor eigenstates: 


'll kqn (j *) 


( *L,n( r A = Jkr ( e -i *' r 

V ei ' l " r U Un( r ) /’ 

(45) 


We have mentioned above that the Bloch states, 
Eq. ( |45| ), are periodic on the interval [0,47r/a) of q- 
points. However, as we will show now, energy dispersion 
£ kqn and wave functions to a q-value in [27r/a, 47r/a) can 
be derived from corresponding quantities in the interval 
[0,27r/a). Therefore, the effective number of spin-spiral 
vectors at which the electronic structure needs to be cal¬ 
culated from first principles is reduced by a factor of two. 

At a given spin-spiral vector q , the Hamiltonian H^) 
has eigenvectors T kqn and eigenvalues £ kqn ■ By sym¬ 
metry, RA) i s identical to the Hamiltonian H [ - q + G ) at 
q + G, where G is a reciprocal lattice vector. Conse¬ 
quently, both Hamiltonians have (i) the same eigenvalue 
spectrum, i.e., £ kqn = £k ' q+Gn and (ii) the same set 
of eigenfunctions such that T kqn = d'fc'q+Gn■ In the 
first principles calculation, these eigenfunctions obey the 


generalized Bloch theorem, Eq. (451, which allows us to 
determine the above crystal momentum k! explicitly: 


’F 


kqn 


(r) = e lkr 


= e 


e* 2 r u I 


l kqn 


kqn 


(r) 

(r) 


ikr 


i— r —i q+G • 

e 2 e 2 


. e 2 


a kqn 
i q + G r ,.-iG r n ,i 


0 ) 


= A k +%)- r 


e 2 e 
». g+ C .r. ~t 


kqn 


( r )y 


. q + G 


U k+§ q+Gn^ 

'^ k+ a q+Gn (r) 


G q+G n 


( r ), 


(46) 


where we defined the lattice periodic functions 


“t ,e^W=“L( r )i 


k+=f q+G n v 

r,4- 


U k+§ q+Gn^ = e 


— iGr 4- 


kqn 


(r). 


(47) 

(48) 


If we change the spin-spiral vector q by G, the Bloch 
state Tfcq„ and its energy £ kqn are moved to a different 
crystal momentum k' = k + G/2. Consequently, Bloch 
states need to be computed only for those spin-spiral vec¬ 
tors which lie in [0,27r/a). 









C. Discussion of the implementation 


D. Computational details 


The evaluation of Mmn\ A), Eq. (34), does not differ 
from the case of standard MLWFs, except that Mmn^ (A) 
needs to be computed for several values of A = q. The 
matrix Mmn\k) in Eq. (351 is given by the overlaps of 
periodic parts Uk qm (r) at neighboring spin-spiral vectors 
q and q + b. If we exploit the generalized Bloch theorem, 
Eq. (45), these overlaps assume the form 


{Ukqm\uk q-\-b n) 

= e±l ^ r (^kqmir))* K{q + b]n( r ) dr ■ 

(49) 


Here, [q + b] is a backfolding of the spin-spiral vector 
q + b to the first BZ, and a =f, J,. The positive (neg¬ 


ative) sign is taken in Eq. (49) for the up component 


(down component) of the Bloch spinor. We describe in 
Appendix |A| the implementation of Eq. (49) within the 
FLAPW method. To reduce the computational burden, 
we can apply the symmetry considerations of Sec. GY! 
to the calculation of the above overlaps. We find that 


M^ G ' b \k)=M^(k+^ 


= lk-j 


and likewise 


(q + G)= M^n ? ,b) (q) = Mmn ' 


(50) 


'(<7), (51) 


where the periodic gauge of the Bloch states in fc-space 
was used. Thus, we can restrict ourselves to the cal¬ 
culation of Mmn\k) and Mmn\q) for spin-spiral vec¬ 
tors in [0, 27r/a). Similarly, the projections (’F kqm\9n) in 
Eq. ( flT| ) need to be computed only for those q which lie 
in this interval. 


Returning to Eq. (35), we have to calculate addition¬ 


ally the auxiliary overlaps. Before, we modify the general 
shape of the auxiliary orbital Cj(£), which was originally 
given by Eq. (14). We choose £ q (£) = e l z^p q (£) such 
that the auxiliary orbital has the same q-period as the 


Bloch states, Eq. (45). The lattice constant in the auxil¬ 


iary space is thus given by a. Then, the auxiliary overlaps 
(pq\p q +b) can be calculated numerically as discussed in 
Sec. |IIIB| As an alternative, we can use the analytic ex¬ 
pression Eq. with r = 27r/(A r q a). 

Projections (' $kqm\gn ) and (( q \h) onto localized trial 
functions q„(r) and /i(£) enter Eq. (37). The FLAPW 
implementation of the former is described in Appendix[A| 
To obtain (C, q \h), we project conveniently onto the single¬ 
well solution, Eq. (151, such that the integral is identical 


to one. However, projections onto different trial functions 
(e.g., Gaussians) can be employed as well. 


Studies of 3d transition metal nanowires indicate that 
spin-orbit effects such as the magnetic anisotropy, which 
we discuss in Sec. |Vj should have a rather small influ¬ 
ence on the electronic structure of Mn atoms due their 
half-filled d shell |28l 125] , Therefore, we neglect spin- 
orbit coupling for the moment. As first step, using the 
one-dimensional version ESS of the density functional the¬ 
ory Jiilich FLAPW code FLEUR 131] . we determine self- 
consistently the electronic density of a one-dimensional 
ferromagnetic linear chain of Mn atoms with a lattice 
constant of a = 5bohr. We employ six local orbitals to 
treat the 3 p core states of Mn. The RPBE parametriza- 
tion of the exchange-correlation potential was used [32] . 
The non-overlapping muffin tin radii and the plane-wave 
cutoff were chosen to be 2.1 bohr and 3.8bohr _1 , respec¬ 
tively. Starting from this charge density, we solve the 
Kohn-Sham equations on a uniform mesh of 8 fc-points 
separately for 16 spin-spiral vectors. 

After that, the information about the wave functions 
at all fc- and q-points is used to compute the necessary 
overlaps and projections. As first-guess trial orbitals g n 
we use three d orbitals and six sp 3 d 2 orbitals for each spin 
direction. The overlaps and projections of the auxiliary 
orbital are derived either numerically or analytically 
as discussed before. A maximal real-space localization of 
the HDWFs is achieved using our extension of the WAN- 
NIER90 code to four space dimensions. Because of the 
metallic character of the magnetic Mn chain, a disentan¬ 
glement [22] °f 18 optimally-connected quasi-Bloch states 
from a manifold of 36 Bloch orbitals is performed. The 
upper bound of the inner, or, frozen energy window is 
2.2 eV above the Fermi energy of the ferromagnetic state 
E F (q = 0) (see Fig.[2|. 

Constructing such HDWFs requires a similar amount 
of computer time as the generation of individual sets of 
MLWFs for all of the 16 spin-spiral vectors. However, we 
emphasize that the single set of HDWFs encodes the com¬ 
plete information of the electronic structure as a function 
of both fc and q in the energy window of interest. 


E. Band structure interpolation 


After constructing HDWFs for the one-dimensional 
chain, we employ these functions in an interpolation of 
the multi-parameter Hamiltonian according to the dis¬ 
cussion of Eq. (38 )-Eq. (41). Figure [2] presents the re¬ 
sults of the generalized Wannier interpolation of the band 
structure compared to the direct calculation. The inset of 
Fig- [2] does not show the usual BZ of crystal momentum 
but a composite BZ combining the crystal momentum 
fc = ke z and the spin-spiral vector q = qe~. According 


to the remark below Eq. (451, the composite BZ is given 


by [—7r/a, 7r/a] x [—27r/a, 27r/a], i.e., it is of rectangular 
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shape. 

A single set of 18 HDWFs allows for an accurate in¬ 
terpolation of the energy bands in the reciprocal (fc, q)- 
space. The path from the T-point (fc = 0 ,q = 0) 
to the X-point (fc = ir/a,q = 0) describes the elec¬ 
tronic band structure of a ferromagnetic Mn chain. The 
standard Wannier interpolation for non-collinear or spin- 
spiral states, which was employed, e.g., in Ref. [33] is al¬ 
ways restricted to high-symmetry lines parallel to T — X 
where q is constant. In contrast, the interpolation based 
on a single set of HDWFs gives access to the electronic 
band structure along any given path in the composite 
BZ. Thereby we can easily compute the electronic band 
structure along the path from the X-point to the point 
M (fc = 7 r/a,q = 27r/a). Along X — M, the crystal mo¬ 
mentum is kept fixed while the texture of the magnetic 
moments changes from the ferromagnetic (q = 0) over to 
the antiferromagnetic state (q = n/a) and back to fer¬ 
romagnetic order {q = 2n/a). Due to the symmetry of 
the band structure discussed in Sec. |IVB[ band energies 
differ at X' and M. The very same set of HDWFs allows 
further for an interpolation of the band structure along 
the diagonal path T — M of the BZ, which is not so easily 
accessible with the standard first principles codes. The 
band energies at T and M are identical due to symmetry 
(cf. Sec. IV B). Overall, the accuracy of the generalized 
interpolation of the band structure is excellent within the 
frozen window. 


F. Real-space visualization of HDWFs 


While MLWFs in magnetically collinear systems with¬ 
out spin-orbit coupling are real-valued [M], they are 
complex-valued functions in non-collinear systems, and 
in the presence of spin-orbit coupling [2T:. In contrast, 
the imaginary part of the HDWFs of spin spirals is neg¬ 
ligibly small such that we can restrict ourselves to a dis¬ 
cussion of the real part of spinor valued HDWFs. 

In the following, we give an argument for the real¬ 
valuedness of HDWFs for spin spirals in absence of spin- 
orbit coupling. We consider the Hamiltonian 

fj 2 

H^ q Ur) = -V 2 + V(r) + B(r — nae z )S n ■ cr 

2m z —' 

n 

_ / rj „—inaq\ 

= ~2- V2 + V ^ + 'E B ( r - na ^{e i ™‘i 0 )’ 

n ' ' 

(52) 


where the first two terms are the kinetic energy and the 
scalar potential. The last term describes the interaction 
with the noncollinear exchange field. The amplitude of 
the exchange field is given by B(r) and its direction is 
given by S n , Eq. (43), within the n-th atomic sphere. It 



Figure 2. Generalized Wannier interpolation of the electronic 
band structure of a one-dimensional Mn chain along high- 
symmetry lines in the composite (fc, q) BZ. The interpolation 
based on HDWFs (red solid lines) is in excellent agreement 
with direct first principles calculations (black circles). Ener¬ 
gies are plotted relative to the Fermi level in the ferromagnetic 
case Ep(q = 0). The thin dotted line indicates the upper 
boundary of the inner energy window of 2.2 eV. 


follows that 





(53) 


If 'F kqn is an eigenfunction of to the real eigenvalue 
£kqn , we can apply a complex conjugation to the corre¬ 
sponding Schrodinger equation and arrive at 


q) (r) (’F kqn (r ))* = £k qn (’F k qn (r))* 


(54) 


where Eq. (53) was used. The complex conjugate of \F k qn 
is an eigenfunction of with energy £k qn - In gen¬ 
eral, eigenfunctions of are labeled by 'S’k-qn such 

that we necessarily need to find (’F fcqn )* = 'F fc /_ q „ for 
some crystal momenta fc and k . From the explicit shape 
of both states dictated by the generalized Bloch theo¬ 


rem, Eq. (45), follows that k' = —fc. We can choose 


the auxiliary orbital cj q of Eq. (141 to obey the relation 
(C q )* = C-q- Thus, the product states in Eq. (18) satisfy 


(<F kqnir , I))* = $-k-qn(r, I) • (55) 

If the unitary matrix satisfies 

(U^)* =U^~ q \ (56) 


the real-valuedness of the HDWFs is implied by Eq. (551: 
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Figure 3. Isosurfaces of a cC^-like HDWF for £o = 0. We 
highlight the icy-plane which is perpendicular to the physical 
chain axis. The up-spin (left) is opposite in sign compared 
to the down-spin component (right), which is of equal magni¬ 
tude. The functions were plotted using the program XCryS- 
Dcn (Ref. [551 ) . 



(x,y) (x,y) 


Figure 4. Isosurfaces of a d xy - like HDWF for z q = 0. The 
auxiliary dimension £ is perpendicular to the highlighted xy- 
plane. Up-spin component (left) and down-spin component 
(right) are of equal magnitude but opposite in sign. 


W'OOri 


N k N X J2 U ™n q) $k q m 

kqm 

2N ^ Nx E (u^ q) *k qm +uL k ’- g) *-k- gm ) 

kqm 

_ 1 _Ft n 

2 ^y N\ ' J \ run ^kqm ' y^mn ^kqmJ J 
kqm 

___» Y" 

AT AT JL / v M mn v kqrn • 
kqm 

(57) 


A very similar argument shows that standard WFs can 
be chosen to be real-valued in some cases: For q = 0 
Eq. (52) describes a magnetically collinear system with¬ 
out spin-orbit coupling, for which Eq. (53) and Eq. ( [55] ) 
simplify into ( H(r ))* = H(r) and (^ fen (r))* = 'F_ fc „(r), 
respectively. The choice Umn^ = ( Umn )* in Eq. |l]) leads 
then to the real-valuedness of the resulting WFs. 

To visualize HDWFs in real space, we first divide a 
given HDWF by its phase at the maximal absolute value. 
Then, one of the four spatial coordinates (:r, y, z or £) 
is kept constant to obtain the three-dimensional plots of 
Fig.[3]and Fig. [4] In general, we find that for a fixed aux¬ 
iliary coordinate £ = £o, the HDWFs (see Fig. [3]) closely 



Figure 5. Scheme of the translation property Eq. (60) of HD¬ 
WFs in the composite lattice (black circles). Both spin com¬ 
ponents of the function Woon are localized in the home unit 
cell ( 0 , 0 ). In contrast, the spin components of Wn 3 n are dis¬ 
placed in the auxiliary dimension with respect to the position 
(ii, 23). Arrows indicate the corresponding distance vectors 

an d 4s,,. ■ 


resemble the first-guess functions of d and sp 3 d 2 charac¬ 
ter visualizing the chemistry of the one-dimensional Mn 
chain. Variations of the value £o do not change the or¬ 
bital character qualitatively. Choosing a constant value 
z = Zq along the chain axis, we present the shape of a d xy - 
like HDWF in Fig. [4] The HDWF extends throughout a 
single unit cell as a function of £ due to the construction 


of the auxiliary orbital, Eq. (14), based on the deep-well 
limit. 

In Fig. [3] and Fig. [4j we show HDWFs in the home 
unit cell, i.e., R = 0 and E = 0. How can we obtain 
the functions Wr=u to finite R and E from those in the 
home unit cell? The spinor components are constructed 
according to 


H kq 


Tg-iq-fS-i) 


XWfcq„(T-)Pq(£) 


(58) 


which follows from Eq. (24), Eq. (451, and the choice 


of the auxiliary orbital. Here, the Bloch-like periodic 
parts U kqn = Y.rn ^mn ] u k qm contain the unitary gauge 
matrix, and <7 =f, 4 .. In the case of the usual WFs, simple 
lattice translations need to be applied to obtain WR n {r) 
to any R, i.e., WR n {r) = W 0n (r — R). Compared to 
the usual WFs, we find from Eq. ( [58] ) a slightly more 
complicated relation between Hoo n (r',£) localized in the 
home unit cell and WRs n (r,£): 

WS 0n (r, £) = WSs n (r, £ + 2E) = W^ 0n (r + R,£± R), 

(59) 

and thus 


WWr, £) = Wg 0n (r - R, £ — 2E T R), (60) 

where the upper (lower) sign is for the up (down) compo¬ 
nent of the spinor. We depict in Fig. [5] the above trans¬ 
lational property of HDWFs for spin spirals. Due to the 
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coupling of fc and 
r (see Eq. (45)), 


q to the same real-space coordinate 
a spin-dependent shift of the spinor 
components occurs for finite R. We can consider the 
distance vector between the centers of Wo 0n (r, £) and 


t a 

l RB.n ~ 




R 
23 ±R 


w^ Sn )-(w; 


OOn 


w, 


OOn 


(61) 


which follows from Eq. (60). If IFoon( r j£) localized in 
the home unit cell at the position (r c ,£ c ), W^ Sn (r,^) is 
centered at (r c +R, g c +2E±R) as shown in Fig. [5] While 
the direct lattice vector R determines the r-center of the 
HDWFs of spin spirals, the center in £-space depends on 
both S and R. 


G. Heisenberg exchange constants and spin 
stiffness 




Starting from the generalized Wannier interpolation of 
the band structure throughout the (k, q)- space, we cal¬ 
culate the dispersion E(q) of the system as the sum of 
occupied eigenvalues for a given value of q. Although 
the energy bands are interpolated nicely using a coarse 
mesh of 8 fc-points and 16 q-points as shown in Fig. [2} 
we find by comparison with direct first principles results 
(see upper panel of Fig. [6]) that an ab initio mesh of 16 
fc-points and 24 q-points is necessary in order to interpo¬ 
late the dispersion 25(g) of the system accurately. This 
means that Bloch functions need to be computed from 
first principles on a (fe, q)-mesh of 16 x 12 if the symme¬ 
try relations from Sec. |IV B| are exploited. 

The value go which minimizes the energy 25(g) de¬ 
fines the ground state of the magnetic system among 
all possible ferromagnetic, antiferromagnetic, and non- 
collinear spin-spiral configurations. Our interpolation 
of E(q) in terms of HDWFs reproduces precisely the 
curve obtained from direct calculation and thereby pre¬ 
dicts the spin-spiral state with go = 0.314 • 2n/a as the 
ground state of the one-dimensional Mn chain at the con¬ 
sidered interatomic distance. The energy difference be¬ 
tween ferromagnetic state and ground state amounts to 
25(0) —25(g 0 ) = 55.2 meV. We further extract the Heisen¬ 
berg exchange constants Jq n by fitting Eq. (44) to the 
HDWF-interpolated dispersion 25(g). The lower panel of 
Fig. [6] reveals that the exchange constants compare ex¬ 
cellently with previous work [25] . 

Wannier interpolation is particularly rewarding in the 
computation of transport properties as crystal momen¬ 
tum derivatives of the Hamiltonian can be taken analyt¬ 
ically [U Him [37]. Such derivatives determine, for ex¬ 
ample, the velocity operator v = V kH^/%. In the case 


Figure 6. Top: Generalized interpolation of 25(g) — E{ 0) of 
the magnetic Mn chain. HDWF-interpolation (solid red line) 
reproduces the direct results (open circles) if the HDWFs are 
constructed using TV*. = 16 fc-points and N q = 24 q-points. 
Bottom: The Heisenberg exchange constants Jon obtained 
by fitting Eq. (441 to the Wannier interpolated energy 25(g) 
(circles) are in excellent agreement with Ref. 1291 (triangles'). 


of spin spirals, the derivative of with respect to q 

can be conveniently obtained from generalized Wannier 
interpolation: 


dH ik '^ 

dq a 


Y,iZ a e' k - R e iq - 3 H(R,E). 

RS 


(62) 


Here, H(R, S) is the matrix of the hoppings H nrn (R, 3), 
and q a and S Q refer to the cr-th components of the vec¬ 
tors q and 3. respectively. Such expressions allow us to 
calculate the second derivative of the energy E(q) conve¬ 
niently, from which we obtain the spin-stiffness A: 


1 d 2 E(q) 
~ 2 dq 2 



(63) 


Necessary details on the implementation of the scheme to 
obtain derivatives of 25(g) are provided in Appendix [E] 
From the evaluation of Eq. (63), we obtain a value of 
A = —174.1 meVxA 2 for the spin stiffness of the one- 
dinrensional magnetic chain in the vicinity of the ferro¬ 
magnetic state. To verify the estimated value for the spin 
stiffness, a polynomial even in q is fitted to the interpo¬ 
lated dispersion near q = 0. We extract a reference value 
of —173.4 meVxA 2 from this fit, which is in very good 
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agreement with the spin stiffness obtained from calculat¬ 
ing directly Eq. (63). 


V. APPLICATION TO THE MAGNETIC 
ANISOTROPY IN A CHAIN OF Mn ATOMS 

A. Introduction 

In this section, we discuss HDWFs for the interpolation 
of the multi-parameter Hamiltonian where rh is 

the ferromagnetic magnetization direction. As an appli¬ 
cation, we consider the magneto-crystalline anisotropy 
energy (MAE), which is the energy difference between 
hard and easy axis of the system. Therefore, we adapt 
our description of the one-dimensional magnetic chain of 
Mn atoms to include spin-orbit coupling. The magnetiza¬ 
tion direction rh = (sin0cos</>, sin0sin<?!>,cos0) is speci¬ 
fied in spherical coordinates by 9 and <j>. Here, we restrict 
ourselves to <j> = 0. Bloch spinors and their periodic parts 
carry a dependence on 0, i.e., '^kdni'r) = e 7kr Uk8n(r). 
We include spin-orbit coupling by the second-variation 
scheme [58] , 


B. Discussion of the implementation 


According to Sec. El C 3[ HDWFs can be constructed 
using projections of the Bloch spinors and overlaps of 
their periodic parts. In the second-variation scheme, 
which we employ to include spin-orbit coupling, the co¬ 
ordinate system in spin space rotates together with rh. 
Consequently, the spinors u k8n (r) = ( u{ gn (r), u l k 0 n (r )) 

and u ke+bn ( r ) = ( u ie+bn( r )i u te+bJ r )) refer t0 Affer¬ 
ent spin-coordinate systems when 6^0. Thus, we need 
to transform the spinors into a common spin-coordinate 
frame when we compute the overlaps of the periodic parts 
at neighboring angles 9 and 9 + b. We use the unitary 
rotation 


X{0) = 



- sin f 
cos | 


(64) 


to obtain all periodic parts UkOn in the same global spin- 
coordinate frame by u k9n = x{0) u k8n- Thar* the over¬ 
laps are given by 

= (^k8m\ U k8+bn) 

= [xHQ)x(Q + b)] a(T , (Uk8m\ u k6+bn ) J 

aa' 

(65) 


where a =t,i- The matrix elements ( u k8m\ u k8+bn) are 


(' u k6m\ u k 


6-\-b n/ 


K[9+b]n( r ) dr , (66) 


where [9 + b\ is a backfolding of the value 9 + b to the one- 
dimensional BZ. These overlaps do not contain an addi¬ 
tional ^-dependent phase as in the cases of the spin spiral, 
Eq. (49), and standard MLWFs. We provide additional 
details and derive corresponding expressions to construct 
the matrix elements within the FLAPW method in Ap¬ 
pendix [d] Apart from the overlaps Mmn\k) we also 
need the overlaps Mmn\ A) with A = 9 (see Eq. (34)). 
However, the calculation of Mmn ^ (9) does not differ from 
the case of standard MLWFs except that several values 
of 9 need to be considered. 


C. Computational details 

We determine the charge density of the ferromagnetic 
chain self-consistently using the computational setup of 
Sec. |IV| but including now spin-orbit coupling in second- 
variation. Based on the converged electronic density, 
we solve on a uniform fc-mesh the Kohn-Sham equa¬ 
tions for each magnetization direction separately. The 
values of 9 are chosen from the range [0, 47 t) as Bloch 
spinors acquire a minus sign upon 360° rotation, i.e., 
ke+ 2 nn = k8n ■ However, symmetry considerations 
analogous to Sec. |IV B| reduce the number of angles 9 for 
which the Bloch functions need to be computed. 

Then, the overlaps Mmn\9) and projec¬ 

tions are calculated. We project onto the same set of 
localized trial functions as in the case of Sec. |IV[ and 
further incorporate the analytical solution for the auxil¬ 
iary orbital. After performing a disentanglement of 18 
optimally-connected quasi-Bloch states from a manifold 
of 36 Bloch states with an inner window up to 2.2 eV 
above the Fermi energy Ep(9 = 0), we generate maxi¬ 
mally localized HDWFs using our extension of the WAN- 
NIER90 program to four dimensions. 


D. Magnetic anisotropy 


The single set of HDWFs is employed to interpolate the 
energy bands £kdn in fe and 9 as described in Eq. ( [38] )- 
Eq (41). In Fig. [7] the energy difference E(9) — E( 0) 
is shown as obtained from such an energy interpolation. 
Compared to the spin-spiral application of Sec. |IV[ the 
HDWFs have to be constructed on a denser ab initio 
mesh of 24 k- and 32 0-points (in [0, 47 t)) to reproduce 
the energy difference accurately. We associate this par¬ 
ticularity with the small MAE of the one-dimensional Mn 
chain of E( 0) — E(rr/2) = 0.217 meV. Thus, a ferromag¬ 
netic magnetization direction perpendicular to the chain 
axis is favored over a parallel orientation as predicted in 
Ref. m 

The uniaxial anisotropy energy can be parametrized 
by E(0) = I\\ sin 2 9, where K\ is the first anisotropy 
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e ( tt ) 


Figure 7. Generalized Wannier interpolation of E(6) — E( 0) of 
the magnetic Mn chain. Irrespective of the small energy scale, 
we can interpolate (solid red line) the energy difference as a 
function of the magnetization direction in very good agree¬ 
ment with the direct calculation (open circles). 


constant, and it follows that [55] 


I<i 


GE{6) 

dO 


0=7I-/4 


(67) 


Generalized Wannier interpolation can be employed con¬ 
veniently to evaluate the above derivative with respect 
to the magnetization direction. At zero temperature, we 
derive in Appendix [E] the expression 


d E{6) 
dO 


y 

Nk 

£ken<E F (9) 


d^kOn 

86 


— y 

N k ^ 


£kGn<EF(0) 


Pk6r 


dH (k ' 6 '> 

dO 


< Pk9r 


( 68 ) 


Here, \(fik6n) are eigenstates of ijWG and the derivative 
of H can be obtained conveniently within the gener¬ 
alized Wannier interpolation scheme: 


80 


yiEe ik R e ie3 H(R, 2 ), 

RE 


(69) 


where H(R , 2) is the matrix of the hoppings H nm {R, 5) 
between HDWFs. The derivative of E(6) at zero tem¬ 
perature, Eq. (68), is the sum of the torques d£ k e n /d0 
which electrons of band n moving through the solid with 
crystal momentum k exert on the magnetization. Us¬ 
ing this approach, we compute an anisotropy constant of 
Ki = —0.224 meV, which agrees nicely with the value 
for E(7 t/2 ) — -E’(O) given above. 


VI. POSSIBLE FURTHER APPLICATIONS 

Above, we have shown that HDWFs can be con¬ 
structed for first principles Hamiltonians and we dis¬ 
cussed applications such as spin stiffness and MAE. In the 
following, we explore within models additional promising 
applications of HDWFs. 


A. Virtual crystal approximation (VCA) 


The electronic structure of non-stoichiometric disor¬ 
dered alloys such as Fe^Coi-a, can be computed within 
VCA where virtual atoms with electronic structure cor¬ 
responding to the concentration x constitute a regular 
lattice [40j. The FLAPW method with a properly ad¬ 
justed number of valence electrons is well-suited to de¬ 
scribe these systems in case of alloys composed out of 
neighbors in the periodic table like Fe and Co ED- % 
computing the electronic structure for several values of 
x , HDWFs can be constructed which describe the multi¬ 
parameter Hamiltonian H^ k ^ x ) for any concentration x. 
Thus, the treatment of alloys such as Fe^Coi-a, on a 
dense mesh of concentrations is simplified. The gauge of 
the alloy Hamiltonians is guaranteed to be smooth due 
to the single set of HDWFs used in the generalized inter¬ 
polation. If one simply mixes the MLWFs for x = 0 and 
x = 1, such a smooth gauge is more difficult to achieve 


Here, we employ a toy model to outline the basic prin¬ 
ciple leaving the implementation into FLAPW for future 
work. We study modulations of the depth of attrac¬ 
tive potential wells at positions Rj = ja which define 
a one-dimensional lattice with lattice constant a along 
the z axis (see Fig. [l] for a sketch of the potential pro¬ 
file). The corresponding one-dimensional single-particle 
Hamiltonian carries a parametric dependence on the vari¬ 
able A: 

H w (z) = ~ (i + asinAjVbX] 0 ^^)’ ( 70 ) 

Rj 


where |a| < 1, and the well function O b R . (z) is defined by 
Eq. (13). The width b and the potential strength Vo are 
chosen such that the three lowest energy bands form an 
isolated group. The Hamiltonian of Eq. ( [70| is diagonal¬ 
ized in a plane-wave basis on a uniform 8x8 (At, A)-mesh. 
The A-points lie in the interval [0, 2n). Necessary over¬ 
laps and projections onto localized Gaussians are con¬ 
structed and the information on the auxiliary orbital 
is derived numerically as discussed in Sec. Ill B Then, 
the WANNIER90 program is used to achieve a maximal 
localization of HDWFs. Figure [8] demonstrates that the 
generalized interpolation reproduces the electronic band 
structure as a function of k and A accurately. 
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Figure 8. Electronic band structure of the multi-parameter 
Hamiltonian Eq. ( |70| ) as obtained by generalized interpola¬ 
tion (solid and dashed lines). The dispersion is depicted as a 
function of the crystal momentum for constant A (left panel), 
and as a function of the parameter A for constant k (right 
panel). The exact results (open circles and diamonds) agree 
nicely with the interpolation for the isolated group of energy 
bands. We have chosen the model parameters a = 3.0 bohr, 
b = 2.9 bohr, a = 0.1, and Vq = 544.0 eV. 


Figure 9. Electronic band structure of the multi-parameter 
Hamiltonian Eq. © as obtained by HDWF-interpolation 
(solid and dashed lines). We present the dispersion as a func¬ 
tion of the crystal momentum at fixed A (left panel), and as 
a function of the parameter A at fixed k (right panel). The 
interpolation is in excellent agreement with the exact results 
(open circles and diamonds). We have chosen a = 3.0 bohr, 
b = 0.5 bohr, and Vo = V/ = 272.0 eV. 


B. Ferroelectric polarization 


In ferroelectrics like the perovskite oxide BaTiC >3 [4B] , 
a relative displacement characterized by the vector A of 
one of the sublattices leads to a change in ferroelectric 
polarization. To determine the value of this change, the 
polarization in the form of MLWF centers or the Berry 
phase has to be computed along a certain path in A- 
space Ml- We can use the HDWFs approach in order 
to interpolate the electronic structure of H along the 
A-path. 

We use the following simple model to describe displace¬ 
ments between sublattices: 


H( x \z) 


d^ 
2to d z 2 


E 

Rj 


v 0 e b R .(z)+v^e b R .(z-T X ) 


(71) 

where t\ = a/2 + 5\ and 5\ = —{b/ 2) sin A describes 
the relative displacement. In addition to the first well 
of depth Vo, which is kept fixed, the same unit cell con¬ 
tains a second well of strength V/ at a variable position. 
Again, we select the well parameters such that the two 
lowest bands form an isolated group. Employing a plane- 
wave basis, we diagonalize the Hamiltonian of Eq. © 
on a mesh of 8 fc-points for each of the 8 parameters A 
chosen uniformly from the range [0,27 t). Having at hand 
the Bloch states T k\n, we construct overlaps and pro¬ 
jections onto Gaussians. After deriving numerically the 


auxiliary orbital as discussed in Sec. |IIIB[ we use the 
WANNIER90 code to establish a maximal localization of 
the HDWFs. The band structure results of the gener¬ 
alized Wannier interpolation presented in Fig. [9] are in 
excellent agreement with exact results. 


For a given value of A the ferroelectric polarization can 
be obtained as sum over the centers of WFs constructed 
from the occupied bands: 


^A = -^ E ( W °n\r\W x n ). (72) 

n(E occ 


Here, e > 0 is the positive electron charge, V is the unit 
cell volume, and \W^ n ) is a standard WF constructed 
according to Eq. ([I]) for H^K The above ferroelectric 
polarization is not unique but defined up to the polariza¬ 
tion quantum. Only polarization changes are unique and 
thus physical. However, to determine unambiguously the 
change of ferroelectric polarization between two points Ai 
and A 2 according to Eq. (721, we usually need to ensure 
a smooth gauge of the Bloch states in A-space. Such a 
gauge is guaranteed if we rewrite the ferroelectric polar¬ 
ization, Eq. (72), in terms of HDWFs. 


How can we now compute P\ within the formalism of 
HDWFs? By performing a Fourier transformation in H, 
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we obtain WFs at A: 


= ^E e “ ifcR ^ X) ^Am(r ,0 


km 


Ca«) 


N k 

= (»«)• 


e ,:fc H ^; A)vI; fcAm(r) 


km 

r\ 


(73) 

Recall that we generate just a single set of HDWFs encod¬ 
ing the electronic structure information in (k, A)-space. 
Consequently, a smooth gauge is automatically built into 
the construction of HDWFs such that also the above 
standard WFs are guaranteed to be smooth in A. In¬ 
serting Eq. ( |73| ) into Eq. ( [72] ) yields 

Px = ~VN~ E 5^e^ (s '- s) (W os „|r|Wo & „) 

Weo COBB' ( 74 ) 

= -t 7 E E elA ' S ^oo„MWo S „), 


nEocc S 

where we exploited 


«>|W 0 A n ) = ^<CaW*MCaW* >, 

-/Va 


(75) 


which follows from (Ca|Ca) = N\. Equation (74) can 
be employed to obtain an interpolated value for the fer¬ 
roelectric polarization at values A that lie between the 
points of the coarse A-mesh used to generate the HD¬ 
WFs. Of course, this works only if the system is insu¬ 
lating along the entire considered A-path. The A-sum of 
P\ evaluates to 

E p * = -f E EE eiA ' H cwwir|Wo S „> 

A nGocc S A 

eN\ 


V 


55 (WoOnMWoors) 


(76) 


which is determined by the centers of HDWFs available 
in the extended WANNIER90 implementation. 

Analogously to derivatives of the multi-parameter 
Hamiltonian discussed in Appendix [EJ we can calculate 
A-derivatives of the ferroelectric polarization, Eq. (74), 
which read 

5Pa -;ee 


•-—< AX-E 


{W 00n \r\W 0 ~ n ). (77) 


d\ a V 

nEocc E 

Here, A a and S a are the a-th components of A and S, 
respectively. Differentiating the ferroelectric polarization 
with respect to the sublattice displacement 6 (which de¬ 
pends on A), we obtain the Born effective charge tensor. 
For the one-dimensional model of Eq. it follows that 


dP s 


dP x d A 


2 dP x 

dS 

6=0 

d\ dS 

A=0 

b d\ 


A=o 


= | E E* S < W OOnMWo~n>. 


(78) 




A 


Figure 10. Treating the tilting angle A as additional vari¬ 
able in the multi-parameter Hamiltonian, we can interpolate 
(solid and dashed lines) accurately the electronic band struc¬ 
ture throughout (fc, A)-space. The dispersion is depicted ei¬ 
ther as a function of the crystal momentum for constant A 
(left panel), or as a function of the parameter A for constant 
k (right panel). Open circles and diamonds indicate the exact 
results. The model parameters are a = 3.0 bolrr, b = 1.0 bohr, 
V 0 = 272.0 eV, and B 0 = 27.2 eV. 


C. Current-induced torques in noncollinear 
magnetic systems 

Current-induced torques on the magnetization (spin 
torques) are thought to play an important role in fu¬ 
ture magnetic memory devices. These spin torques result 
from the exchange of angular momentum between two 
magnets of distinct magnetization direction (spin trans¬ 
fer torques) j44H46] . or between spin and lattice (spin- 
orbit torques) [lOj |TT] l47ti49j . The spin-orbit torques 
can depend strongly on the magnetization direction m 
We expect that HDWFs provide a convenient scheme to 
extract this directional dependence. 

To demonstrate that the generalized interpolation of 
the Hamiltonian with respect to isolated spin moment 
rotations in real space can be tackled with the formalism 


of HDWFs, we modify Eq. (70) to describe two magnetic 


“atoms” separated by half the lattice constant s = a/2. 
And while we keep the orientation of one of the atoms 
fixed, i.e., h\ = e x , the direction of the other moment, 
h 2 = (cos A, sin A, 0), is tilted by the angle A. The result¬ 
ing single-particle Hamiltonian assumes the form: 


hW(z) =- L ^ - v ° £ K <*>+•a - »> 

" Ri 


+ Bq 55 [©Rj {z)hi + <d b R . ( 


z - s)n 2 


bV 


( 79 ) 
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this spin part: 


/W(0\ 


E 

A 

E 





/ e iA (! + s\ 


(80) 


where the phases e 1 ^ guarantee the orthogonality, and 
the gauge freedom is represented by cj>\. It follows that 
W^(£ — 2) = —W^(£). Consequently, the components of 
the spinor are opposite in sign and additionally shifted 
by two lattice constants in £-space. 


D. Mixed Berry curvature 


Figure 11. Real-space distribution of the up (left) and down 
(right) component of the HDWF associated with spin mo¬ 
ment rotations. Red balls refer to those moments ho which 
rotate with A, and blue balls represent the positions of the 
fixed moments fi i. Black lines indicate contours of constant 
function value. The functions were plotted using the program 
XCrySDen (Ref.®. 


where Bq is the strength of the exchange potential, and <x 
is the vector of Pauli matrices. Diagonalizing the matrix 
Eq. (79) in a plane-wave basis on a coarse 8x16 (fc, A)- 
mesh allows us to extract the Bloch spinors, and to calcu¬ 
late the matrices necessary to apply the WANNIER90 min¬ 
imization to HDWFs. Recalling that the Bloch spinors 
acquire a Berry phase of ir upon rotating by 360°, we have 
to choose the A-points uniformly in the interval [0, 47 t) 
(see also the remark below Eq. ([45])) . The auxiliary or¬ 
bital is taken as Ca(£) = 

As shown in Fig. [To] HDWFs succeed in the precise in¬ 
terpolation of the band structure of the family of Hamil¬ 
tonians Eq. (79) throughout the composite BZ of fc and 
A. For this model, we also present in Fig. 11 the lo¬ 
calized real-space distribution of one of the low-energy 
HDWFs. Remarkably, the spinor-valued HDWF turns 
out to be purely real. The individual components are 
both centered in the potential wells at s = a/ 2 , where 
the exchange field is rotated. However, the down com¬ 
ponent is displaced by +2 lattice constants along the £ 
direction compared to the up component. 

In the following, we give a simple argument for the shift 
in the coordinate £ between up and down components of 
the HDWF shown in Fig. El If we consider the deep- 
well limit Vo —► oo, the two atoms in the unit cell do not 
hybridize and decouple completely. Thus, the problem 
is equivalent to finding the lowest-energy solution for a 
chain where all moments rotate with A. The spin part 
( e *A/ 2 ^ _ e -» a/ 2 ) gyoP a solution describes the rotation 
around the z axis of a spin pointing initially in the — x 
direction. We can perform a Fourier transformation of 


Recently, the mixed Berry curvature in (fc, m)-space 
has been found to be important for spin-orbit torques, for 
the Dzyaloshinskii-Moriya interaction and for the charge 
of skyrmions mmmm- This mixed Berry curvature 
is given by 


( / <9u gl 

D" (fc, m) = - 2 e* • m x 3 7 kmn 


dm. 


du sl 

UUj kmn 


dk.j 


(81) 


where i and j are Cartesian directions and e, is the unit 
vector in the i-th Cartesian direction. If an electric field 
E is applied to a ferromagnet with broken inversion sym¬ 
metry, the torque 

r = - ] ^^^eO”.(fe,m)e i F J -, (82) 

k kn ij 


acts on the magnetization due to the mixed Berry cur¬ 
vature D™(fc,rh) EHini. Using spherical coordinates 
to express the magnetization direction such that rh = 


(sin d cos</>, sin 0 sin ij), cos 6), we can rewrite Eq. (81) as 


n«(fc,M) = - 2 e* 


f . E du te4,n 
e' 


86 


8v Bl 

UUj k8<t>n 


dkj 


1 


—— 
sin 6 \ o<p 


/ du ie4>n du te<t,n ' 


dkj 


(83) 


If we construct HDWFs for the Hamiltonian we 

can use the generalized Wannier interpolation in order 
to evaluate Eq. (83). In Sec. [V] we demonstrated that 
HDWFs can be constructed for . It is straight¬ 

forward to extend the scheme of Sec. m to allow for 
variation of both 6 and <j>. How to obtain the deriva¬ 
tive | dukn/dkj) from the standard MLWF interpolation 
is discussed in detail in Ref. |2 : . The above derivatives 

I du tefJ dk j)’ \ du l E/- 961 )’ and \ du lwJW) are emu¬ 
lated in the HDWF-interpolation scheme in a similar way. 
We suppress the superscript “gl” in the following. 
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The periodic parts of the Bloch-like functions are ob¬ 
tained from the HDWFs by Fourier transformation: 


\<P, 


kd(pn ) bt 

-EE 


e -»*»■(»— R) e -i8(£o-£e) 


R 


(84) 


x e 


^0(^0 1=1 4>) 


IkF, 


Q iHj (f) Ti 


>• 


In order to acquire the periodic parts of the eigenfunc¬ 
tions of (cf. Eq. (40)), we need to apply an 

additional unitary matrix yww (cf. Eq. ©): 


l^kGtfin) ~ E ^fikec/inJVmn 9 ’^ — I U kBc\>nPs^ > 


where |u fee<M ) = E m Accordingly, the 

0-derivative of the latter is given by 



where the derivative of the unitary transformation is 
written as dV^’^/dB = y(M,0)£)(M,«. Introduc¬ 
ing the abbreviations A = (8,<f>), £ = (f;@. Ctf>), and 
S = (Eg, E^), we can use for the first term that 


( 


U ke<t>n 


^ u kG(j>m \ ST 'ST „ik-(R-R’)„i\-( S-l 

de ~2^ 2^ e 

' PS R'.w.' 


) 


x (FFr's'„| [—?'(£e — Sg)] |Wjiam) 

R3 H'H' 

(87) 


which follows from Eq. ( [84] ) and (p a |Va/ 9 a ) = 0- The 
matrix elements {Wr"=’ n \^,e\^RSm) can be computed by 
generalizing Eq. (27): 


(wwjsiww) = 


N k N } 


■E 

kX 


^-ik-(R-R') -iA-(H-S') 


X (^fcAnl^Al^fcAm) • 


( 88 ) 


Similar off-diagonal matrix elements are available in the 
WANNIER90 code for the case of standard MLWFs. They 
are obtained by approximating the gradient by finite dif¬ 
ferences [22] : 


(W 0 n\r\W Rm ) 


i 

Nk 


kb 




(89) 

It is straightforward to generalize 
case based on the overlaps in Eq. 


Eq. (|89J) for the HDWF 
([Ml) and Eq. (l35|). For 


the second term in Eq. (86) we can use [2] 


lj(k,9,<l>) 

ran 


/ 


\ 

\rkO(f)m 

dG 

rkG(f>n j 


&k6(l>n £kdc/>m 


if n ^ m 


I 0 if n = to 

(90) 

which can be evaluated for any (fc, 0, <j >) from general¬ 
ized Wannier interpolation. Analogously, the derivatives 
I d u ke<t>n/®kj) and |<9u fce 0 ra /<9<(>) are constructed within 
the formalism of HDWFs. 


VII. SUMMARY 


We introduce the concept and formalism of higher¬ 
dimensional Wannier functions (HDWFs) to describe 
the electronic structure of multi-parameter Hamiltonians 

if( fc A) ; w h ere \ 

is an external periodic parameter. The 
introduction of an auxiliary space £ solves the fundamen¬ 
tal problem of non-orthogonality of usual Bloch states 
in such a situation. Analogously to maximally-localized 
Wannier functions, we define HDWFs as Fourier trans¬ 
formations of higher dimensional product states carrying 
a dependence on k and A. A minimal and accurate inter¬ 
polation of multi-parameter Hamiltonians is established 
using HDWFs. The implementation of the necessary ma¬ 
chinery for the construction of HDWFs from ab initio 
within the FLAPW method is discussed. In order to 
achieve a maximal localization in the extended space of 
r and we adapt the WANNIER90 program. 

The application of the formalism to a one-dimensional 
Mn chain with the spin-spiral vector as an external pa¬ 
rameter reveals an excellent agreement with direct first 
principles calculations, and enables the simplified extrac¬ 
tion of Heisenberg exchange constants and spin stiffness. 
Treating the direction of the ferromagnetic magnetiza¬ 
tion in real space as an external parameter, we are able 
to apply the HDWFs machinery to compute the magneto¬ 
crystalline anisotropy energy of the Mn chain. Although 
the corresponding energy scale is very small and thus 
more difficult to capture, HDWFs interpolate accurately 
the energy E(6) as a function of the magnetization di¬ 
rection. We outline various physical problems to which 
HDWFs could be applied efficiently, e.g., disorder treated 
within VCA. We emphasize further the advantages asso¬ 
ciated with the evaluation of linear response coefficients 
such as AHE and spin torques. A formula for the gener¬ 
alized interpolation of the ferroelectric polarization along 
any insulating A-path is provided. Finally, HDWFs could 
prove useful in the topological characterization of com¬ 
plex multi-parameter systems as they allow for the sim¬ 
plified evaluation of mixed Berry curvatures. 
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Appendix A: Details on the FLAPW implementation 


(k b) 

Expressions for the overlaps M„ lr { between periodic 
parts of the Bloch states at neighboring crystal momenta 
and projections iSi were already derived for an imple¬ 
mentation of MLWFs within FLAPW |27j . The evalua¬ 
tion of Eq. (35) requires additionally the construction of 
the overlaps 


M^\k)=Y J (ul qm \ul q+bn ) (Al) 

a 


with the vector b = bb connecting the two q-points, and 
a =t,4~ Because of the real-space partition into muffin 
tin spheres (MT) and the interstitial region (INT), these 
matrix elements decompose: 


M^\k) = M^ b \k) + £ M^\k) 

INT z —' 


MT U 


(A2) 

Here, /r labels the different atoms in the unit cell. Fur¬ 
ther contributions arise in film calculations (cf. Ap¬ 
pendix [B]), the study of one-dimensional geometries (cf. 
Appendix 0, and when the FLAPW basis set is sup¬ 
plemented with local orbitals. Due to the generalized 
Bloch theorem, Eq. (451, the overlaps in Eq. ( |A1[ ) can be 
rewritten: 


M± b \k) = J2 [ e±l| ' r (n qm (r)Yn iq+b] u(r) dr, 

(A3) 

where the upper (lower) sign is associated with the up- 
spin (down-spin) of the Bloch states. The expression [q] 
refers to a backfolding of the momentum q into the first 
BZ by the subtraction of a reciprocal lattice vector G(q), 
namely, [q] = q — G(q). As a consequence of the doubled 
BZ of q-points, G{q) is twice as large as a usual reciprocal 
lattice vector (see also remark below Eq. ©). 

Within the muffin tin sphere centered around the /r-th 
atom, which is located at the position plane-waves 
do not succeed in describing the physics in presence of 
the singular atomic potential. Thus, the Bloch states 
are expanded in terms of radial solutions ui of the scalar 
relativistic equation at band-averaged energies, related 
derivatives with respect to energy tiq, and the spherical 
harmonics Yj, where L = (l,l z ) represents the set of an¬ 
gular momentum quantum numbers. Accordingly, the 


(A4) 

+Kn ( fc > y L(r„) ■ 

Here, a^ and are expansion coefficients in the /i- 
th muffin tin and the position relative to the nucleus is 
denoted as = r — r M . If we employ the Rayleigh 
expansion 


e^ b r = 4ne^ ib -^ ^T(t1)^M)Yl(&) (Y l (»>))* 

(A5) 

of the plane-wave factor in Eq. (A3) into spherical har¬ 


monics, the muffin tin contribution to the overlaps be¬ 
tween periodic parts assumes the form 


M^ b \k) 


MT„ 


= 4 

a 

x <?))* a%T n (k, [q + b])t%[ LL {b, a) 

LL' 

+ ?))* [<? + 

+ <?)f a L'n( fe ’ [<1 + b ]Ki LL ( 6 ’ cr ) 

+ q ))* K’>l( k Aq + b ))t22 LL ( b w) • 


(A6) 


Here, the radial solutions, their energy derivatives, and 
the spherical Bessel functions ji enter through the t- 
coefficients defined as 


t‘if L {b 1 a) = Y J QLL'L"{ b ) 

L' 

X / r ^ 1 ' ("2”) u ^ ,,T ^) u c <T ( r ' Lt ) dr < u ’ 

ti2 L L ( b iCr) = ^Q LL ' L "{b) 

L' 

X J (jy) K'‘ T ( r ») u i"' T ( r ») dr », 


(A7) 


(AS) 


and likewise for < 2 i and ^ 22 - If we choose a uniform 
Monkhorst-Pack grid to sample the BZ of spin-spiral pa¬ 
rameters, the above integrals become independent of the 
q-point such that they may be calculated once and for 
all at the very beginning. The abbreviation 

Qlul" i b ) = i l (±1)* Y L '{b)GLL’L" (A9) 


incorporates the Gaunt coefficients Gll'L", which are 
given by 

G ll , l »= J Y L (r ll )(Y L ,{r lt ))*(Y L ,,(r li ))' dH. (A10) 

The expressions above are easily extended when local or¬ 
bitals are employed in the basis set. 
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In FLAPW, the Bloch states are expanded using plane- 
waves with reciprocal lattice vectors G in the interstitial 
region. Thus, the wave function assumes a form in line 
with the generalized Bloch theorem: 

n qn {r) = ^ E c g(^ 9, n)e< k ^ +G ) r . (All) 

G 


Defining the Fourier transformation of the step function 
©int cutting out the interstitial region by 


1 


®G=y I e 
V J INT 


1 


d r=^J e~* Gr ©int (r) dr , 

(A12) 

we can write the interstitial contribution to the overlap 
elements of the periodic parts at neighboring spin-spiral 


parameters, Eq. (A3), as 


GGV 


■Mmn\ k ) TM = E (cG( fe i9,™))*CG'(fc, IQ + b \ ,™) 
int 

(A13) 


X ©^G(q + b) 
+ 2 


+G-G' ’ 


The shapes of the overlaps Eq. (A6) and Eq. (A13) 
differ slightly from those of the Mmn ■* contributions de¬ 
scribed in Ref. ED First, the expansion coefficients carry 
a new dependence on the spin-spiral vector q. An ad¬ 
ditional spin-dependent sign arises from the generalized 
Bloch theorem in Eqs. (A6), (A9), and (A13). Finally, 
the vectors b and G(q+b) occur both with a factor of 1/2 
in Eqs. (A6), (A13), and the definition of the f-integrals. 

To construct first-guess HDWFs, the projections of the 
Bloch states onto localized trial orbitals g n have to be 
evaluated within FLAPW according to Eq. (37). These 
trial orbitals are chosen to be zero everywhere except for 
the /it-th muffin tin sphere to which the corresponding 
first-guess should be associated. The expansion coeffi¬ 
cients in g n {r) = J2L c LnUi{r /J ,)Y L (r IJ ,) control the angu¬ 
lar character of the trial functions m- The radial func¬ 
tion ui can be chosen, for example, as the first principles 
solution Ui to the radial Schrodinger equation. Then, 
projections {^kqm\gn) are computed according to 


E 

La 


K^(k,q))*cl n 


( & Lm( fc ’9)f C L 


r l ^Mu/Mdrv 
r l ^{r^uiir^dr^ 


(A14) 


if the orthogonality of the spherical harmonics is ex¬ 
ploited. Except for the q-dependence of the expansion co¬ 
efficients, these expressions are similar to those described 
in Ref. 1271 for standard MLWFs. 


Appendix B: Vacuum contribution to the overlaps 
A imn^k) in film calculations 

In the study of two-dimensional geometries using the 
film implementation of the FLEUR program, an additional 


contribution to the matrix elements in Eq. (A3) occurs as 


a consequence of the presence of two semi-infinite vacua 
& The Bloch states in each of the vacua, which extend 
from —oo to —T >/2 as well as Z?/2 to oo, are represented 
by 

^, 9|l «W=EC( fe ||,q||,z)e< h » T ^ +G »K (Bl) 


Here, fcy and qy are both considered to lie within an 
according two-dimensional BZ associated with the film 
plane, which is supposed to be perpendicular to the z- 
axis. The function 


( fc ll > 9|| > *) = a G||n( fc || . 9)1 )«Gh ( fe !l ’ 9j| , z) 
+ ^G || n( fe |h9||)MG || ( fc |b9||,^) 


(B2) 


includes the one-dimensional solutions of the Schrodinger 
equation in the corresponding vacuum region uc n and 
their energy derivatives ■ For convenience, the ab¬ 
breviations 


£g"g' ( fc ll . 9 1 |, [9 1 | + b],z) = (v£f( fc ll > 9 1 | , zj) 


x ^G' ff ( fc II > 9y + b ,z) 


(B3) 


and = G||—G|| =FGy (q||+6)/2 are introduced. Conse¬ 
quently, the contribution of the vacuum extending from 
T >/2 to 00 to the overlap matrix elements between peri¬ 
odic parts, Eq. (|A3|), evaluates to 


Mmn’ b \k\\) 


FILM 


AE sii%, 


r°° ± f g »t’ a ii+»L 


Giig:, 


I'D/2 


x$%&Ak hqp [q n +li\,z)dz : 


(B4) 


where the unit cell area with respect to the film plane 
is denoted as S y. The other contribution from the sec¬ 
ond vacuum region is derived analogously. Compared to 
the contribution to the usual overlaps Mmn \ Ref. ED 


the function i/X of Eq. (B2) carries a dependence on 


< 7 ||. Additionally, the reciprocal lattice vector G(qn + b) 
occurs with a spin-dependent sign and a factor of 1/2 in 
the definition of Q\\ and Eq. (B4|. 


Appendix C: Vacuum contribution to the overlaps 
Mrnn\k) in one-dimensional calculations 

The density functional theory code FLEUR treats one- 
dimensional systems as cylinders with radius R vac em¬ 
bedded in surrounding vacuum [ 50] , The cylinder axis 
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points along the z direction. Using cylindrical coordi¬ 
nates in real space r = (z , r, <j >) and reciprocal space 
G = (G z ,G r ,Gt), we express the single-particle wave 
function in the vacuum as 

; (C1) 
p 

where as well as q z are drawn from a one-dimensional 
BZ, and the integer p labels the cylindrical angular har¬ 
monics. The variable P denotes the set of ( G z ,p ) with 
respect to which the summation is performed. Radial 
solutions Up to the Schrodinger equation in the vacuum 
region and related energy derivatives up enter the ex¬ 
pression through 

^p a {kz,qz,r) =a^ a {kz,q z )up(k z ,q z ,r) 

+bp a {k z ,q z )up(k z ,q z ,r). 

For convenience, the abbreviations 


(C2) 


(C3) 


Pppf(k z , qz, [q z + b\,r)= °(k z , q z , r))* 
xippf (k z ,[q z + b\ ,r) 

and Q z = G z — G' z ^fG z (q z +b)/2 are introduced such that 
the corresponding overlap elements, Eq. (A3), associated 


with the presence of the vacuum assume the form 


0D =EE 


[«, + &], r) 


pp, 


xe -i5*z e ±i 


Jpp, 

VAC 

G|| (<2z+&) 


e *(p'-p)0 dr. 


(C4) 

Here, r || shall refer to the x- and y-conrponent of the real- 
space vector r = (r\\,z) and similar for the reciprocal 
lattice vector G\\{q z + b), which shifts the momentum 
back into the first BZ. Exploiting then the plane-wave 
expansion into cylindrical coordinates 

e =FiG - r = e TzG * 2 ^ i p (Tl) p e Tip{ ^ G )j p (G r r ), (C5) 

p 

we arrive finally at the vacuum contribution to the over¬ 


laps of periodic parts at neighboring q, Eq. (A31, in case 
of one-dimensional calculations: 

p'-p ap-p' p i(p'-P)<t>a(qz+b) 


^™ ) fe) OD =EE(Ti) p ' _p * : 


pp/ 


x£Sg z 


x/3 


P p, 


I 'f'Jp'—p 

/i?v ac 

’(k z ,q z , [q z - 


G r (q z + b)r 
2 

b }, r) dr . 


(C6) 

Here, J p represents the cylindrical Bessel function of or¬ 
der p, and £ = 2nT with the lattice constant T along 
the axis of translational invariance. In contrast to the 
implementation of the usual overlaps Mmn \ Ref. HI 
a spin-dependent sign arises from the generalized Bloch 
theorem in Eq. (C6). The lattice vector G(q z + b) occurs 


further with an additional factor 1/2 in the argument of 
the cylindrical Bessel function, and the definition of Q z . 


Appendix D: Calculation of Mmn\k) within FLAPW 


Knowledge of the overlaps between periodic parts of 
the Bloch states at neighboring angles 9 and 9 + b is 
required to construct HDWFs when magnetization di¬ 
rection plays the role of the additional external parame¬ 
ter. Within the second-variation scheme [55] used in this 
work, the spin quantization axis of the wave functions is 
identical to the magnetization direction, which we char¬ 
acterize by an angle 9. Using the rotation 


x(Q) = 


(Dl) 


, sm: 


COS ; 


we transform therefore all wave functions to the very 
same global frame in order to evaluate the overlaps 

AfE'W = EKilKIh.) 

a 

= + b )\aa’ ( U k6m\ u k8+bn) ■ 

aa' 

(D2) 

Here, the periodic part ui-g n in the local coordinate frame 
was transformed to the global one b y = x(®) u k 9 n, 
and a =t,4~ Keeping in mind Eq. (D2|, we present in 


the following the necessary FLAPW expressions for the 
calculation of the overlaps in the local spin frame. 

The standard expansion of the wave function into plane 
waves is used in the interstitial region with the expansion 
coefficients carrying now a dependence on the angle 9: 




( r ) = ^E c G(M,ny< fc+G) U 

G 


kOn 


(D3) 


Thus, the overlaps of lattice periodic parts in the local 


coordinate frame, Eq. ( 661 , assume the form 


kOm \^k 6+b n 


INT 


GG' 


{c a G (k, 9 , to))* Cq,( k, [9 + b ], n)0 G _ 


G’ 


(D4) 


where 0 G has been defined in Eq. (A12). Compared to 
the implementation of the usual overlaps Mmn ' 1 , Ref. EH 
only reciprocal lattice vectors G and G 1 enter 0 G above. 
Thus, we can arrive at the shape of the above overlaps 
by formally setting G(q + b ) to zero in Eq. (A13). 

In contrast to Eq. (|A4|), the coefficients of the expan¬ 


sion of the muffin tin wave functions depend on 9. Ac¬ 
cordingly, the Bloch state in the local spin-coordinate 
frame is given as 


^ k9n( r ) 


MT„ 


L 


[ a Ln( k ^) U l ,,T ( r p) 


(D5) 


+ b Ln( k ’ 9 ) il l’ <T ( r p)] Y L{rp.) 
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where L stands for the set of angular momentum quan¬ 
tum numbers (l,l z ). The overlaps between the lattice 


periodic parts, Eq. ( 661 , are evaluated using the orthog¬ 


onality of the spherical harmonics to yield 


MT„ 


i u ke m K Q+bn)\ 

= £ [( a Lm( fe ’ °)T a Ln ( fe ’ \ e + b }) t li L (<7, <*') 
L 

+ (fc, [0 + b])t^ L (a,a') 


d )T a Ln ( fc > i 6 + CT ') 

d )T Kn \° + & ]) i 22 L (^ *') 


(D6) 


where the coefficients ty represent integrals of the radial 
solutions and their energy derivatives: 

tii L (v, <?') = J r l u i’ a ( r fi) u i ,Cr Ov) dr V - ( D7 ) 

*12 L {<T> a ') = j r l U l’ rT ( r ^ il l’ <T M dr M - ( D8 ) 

and likewise for f 2 i and t 22 . Compared to App endix |A| or 
the implementation of the usual overlaps Mmn ' 1 , Ref. [27, 
the above t-integrals are simplified as they do not con¬ 
tain the Gaunt coefficients. Formally, we can obtain, for 
example, a) from Eq. (A7) when setting b to zero. 

If we consider the application to the one-dimensional 
magnetic chain discussed in the main text, an additional 
contribution arises due to the presence of the vacuum 
(cf. Appendix [Cj. The wave function is expanded in the 
vacuum region as 


n z9n (r) = Y / r P ’ ,T (k z ,e,r)e i ^e i ^ 


~\~G z )z 


(D9) 


with P = ( G z , p) representing the set of the integer p and 
the plane-wave vector G z , and further 


ipp a (k z , 0 , r ) =a“(t, 9)up(k z ,6 , r) 

+bp a (k z ,6)up(k z ,d,r). 


(DIO) 


Here, up and up refer to the radial solutions of the 
Schrodinger equation in the vacuum region and their en¬ 
ergy derivatives, respectively. Consequently, the vacuum 
contribution to overlaps of the periodic parts in the local 
frame, Eq. ( |66| ), assumes the form 

(^kdm\^k 6-\-b n) | OD — 

/»oo 

= V’p' 7 (kz,[0 + b],r)dr, 

p Jflvac 

(Dll) 


where £ = 2ttT with T as lattice constant measured 
along the z direction, and R vac is the radius of the one¬ 
dimensional geometry under consideration. Unlike the 


case of the usual overlaps Mm'n ^, Ref. na no cylindrical 
Bessel function occurs in the above radial integrals. The 
formal shape of such overlaps can therefore be obtained 
by considering G{q z + b) = 0 in Eq. (C6). 


Appendix E: Derivatives of the multi-parameter 
Hamiltonian with respect to the additional 
parameter A 


The Wannier interpolation scheme provides an ele¬ 
gant means of performing analytically crystal momentum 
derivatives of the Hamiltonian, which enter the calcula¬ 
tion of properties such as the AHE or other transport 
coefficients d is na mzi- We are able to compute anal¬ 
ogously derivatives of the multi-parameter Hamiltonian 
jf{k,\) re spect to an additional external parame¬ 


ter A, starting from Eq. (40) of the generalized Wannier 
interpolation: 

dH 
d\ 0 


= £* 


rr t rk-R'iX-’E. 


H(R, S), 


(El) 


PH 


and 


d 2 H 

dX a dXp 


= £ 


= Ppe ik R e lX S H{R, 3), (E2) 


RS 


where H(R, S) is the matrix of the hopping elements 
H nm {R , 3) between HDWFs, and Aq, and E a refer to the 
a-th components of the vectors A and S, respectively. To 
simplify notation, we suppress the explicit dependence 
of on k and A here and in the following. The 

above equations may be particularly fruitful in accessing 
accurately Berry connections and curvatures. 

We employ such expressions to determine the first and 
second derivatives of the energy E{ A) with respect to the 
external parameter A. Based on the Fermi-Dirac distri¬ 
bution function f(y) with y = Ep(\) — £kXn , the energy 
of the system is defined by 

EW = ^Y, £ ^f(y), (E3) 


with the Fermi energy Ep{ A), and it follows that 

d a E( A) = [d a £k\nf(y) + £k\nd a f(y)\ , (E4) 

k kn 


and 

d a dpE( A) = -jr- 'Yyd a d/ j £k\nf{y) + £kXnd a dpf(y) 

kn 

~\~da.£kqn^jH f {y) T dfi£kqndaf X?/)] ; 

(E5) 


where the notation d a = d/dX a was introduced. We can 
obtain the derivatives of the band energies, which enter 
these equations, by using Eq. (El) and Eq. (E2): 


$a£k\n — {^PkXn\^aR\^Pk\n) j 


(E6) 
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and 


dcn.dp£'k\n — {^k\n\dad(3H\ ( ftk\n) 

(^k\n\daH\(pk\rri} (^feAm| &fiH\(fk\n) 


+ 2 ^E 

m^n 


£kXn ~ £k\m 


(E7) 


where the second contribution can be derived from first 
order pertubation theory. The states \ipk\n) are the 
eigenvectors of the multi-parameter Hamiltonian H^ k ’ X \ 
Evaluating d a f and d a d p f in Eq. (E4) and Eq. (E5) re¬ 
quires knowledge of the derivatives of the Fermi energy 
Ep{ A). To obtain analytically the necessary informa¬ 
tion, we invoke the total number of electrons in the sys¬ 
tem, N( A) = N ^ 1 Sfcn f(y)> which is a constant, i.e., 
d a N(X) = 0. First derivatives of the Fermi energy are 
accordingly given by 


d a E F ( A) = 


E 

_ kn 


df(y) 

dy 


-1 


E , (E8) 


kn 


dy 


where the term J2kn df{y)/dy is a measure for the den¬ 
sity of states at the Fermi level. The second derivatives 
of the Fermi energy assume the form 


d a dpE F { A) = 
d 2 f{y ) 


E 

_ kn 


df(y) 

dy 


-1 


E 

kn 


df(y) 

dy 


d a dp£, 


k\n 


dy 2 


(■ d a Ep(X ) — d a £k\n){dpEp(X) — d p £ k xn) 


(E9) 


which is easily found from the condition d a dpN( A) = 0. 

At zero temperature, Eq. (E4) and Eq. ( |E5| ) simplify. 
From the condition d a N(X ) = 0 follows that 


d a E(X) = — E f(y)d a £kxn 


kn 


— jy ^ ^ Q(l/) {^kXnldgHl^pkXn) ? 


(E10) 


kn 


with Heaviside step function @(y). Likewise, we obtain 
d a dpE(X) = 2- E (f(y)d a d 0 £kxn + d a f(y)d p £ k xn) 


kn 


kn 


— jy ^ ^ Q(y) | {Pk\n\dctdf3H\(pk\ n ) 


m^n 


{(Pk\n\&aH\(Pk\ni) (^kXm 1^/3-^l^fcAn) 


1 

Nk 


£kXn ~ £kXm 

^ ' fi(y) (d a Ep( A) d a £kXn ) ( t PkXn\d p H\ipkXn ) - 

kn 

(Ell) 

To calculate accurately the derivatives of the energy 
E(X) given by Eq. (E4) and Eq. (E5), we implement 


the above scheme based on the hoppings. We are able 
to derive from generalized Wannier interpolation basic 
properties of the system, for example, the spin stiffness 
or the anisotropy constant. 


* j. hanke@ fz-j uelich.de 

[1] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and 
D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012). 

[2] X. Wang, J. R.. Yates, I. Souza, and D. Vanderbilt, Phys. 
Rev. B 74, 195118 (2006). 

[3] J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, Phys. 
Rev. B 75, 195121 (2007). 

[4] Y. Yao, L. Kleinman, A. MacDonald, J. Sinova, T. Jung- 
wirth, D.-s. Wang, E. Wang, and Q. Niu, Phys. Rev. 
Lett. 92, 037204 (2004). 

[5] R. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 
(1993). 

[6] D. Vanderbilt and R. King-Smith, Phys. Rev. B 48, 4442 
(1993). 

[7] R. Resta, Rev. Mod. Phys. 66, 899 (1994). 

[8] M. Lezaic, P. Mavropoulos, G. Bihlmayer, and S. Blrigel, 
Phys. Rev. B 88, 134403 (2013). 

[9] F. Freimuth, S. Bliigel, and Y. Mokrousov, J. Phys.: 
Condens. Matter 26, 104202 (2014). 

[10] F. Freimuth, S. Bliigel, and Y. Mokrousov, Phys. Rev. 
B 90, 174423 (2014). 

[11] K. Garello, I. M. Miron, C. O. Avci, F. Freimuth, 

Y. Mokrousov, S. Bliigel, S. Auffret, O. Boulle, 
G. Gaudin, and P. Gambardella, Nature Nanotech. 8, 
587 (2013). 

[12] H. Kurebayashi, J. Sinova, D. Fang, A. C. Irvine, T. D. 
Skinner, J. Wunderlich, V. Novak, R. P. Campion, B. L. 
Gallagher, E. K. Vehstedt, L. P. Zarbo, K. Vyborny, A. J. 
Ferguson, and T. Jungwirth, Nature Nanotech. 9, 211 
(2014). 

[13] E. Roman, Y. Mokrousov, and I. Souza, Phys. Rev. Lett. 
103, 097203 (2009). 

[14] F. Freimuth, R. Bamler, Y. Mokrousov, and A. Rosch, 
Phys. Rev. B 88, 214409 (2013). 

[15] P. L. Silvestrelli, N. Marzari, D. Vanderbilt, and M. Par- 
rinello, Solid State Commun. 107, 7 (1998). 

[16] M. Posternak, A. Baldereschi, S. Massidda, and 
N. Marzari, Phys. Rev. B 65, 184422 (2002). 

[17] Y.-S. Lee, M. B. Nardelli, and N. Marzari, Phys. Rev. 
Lett. 95, 076804 (2005). 

[18] H. Abu-Farsakh and A. Qteish, Phys. Rev. B 75, 085201 
(2007). 

[19] V. Anisimov, D. Kondakov, A. Kozhevnikov, I. Nekrasov, 

Z. Pchelkina, J. Allen, S.-K. Mo, H.-D. Kim, P. Metcalf, 
S. Suga, et al. , Phys. Rev. B 71, 125119 (2005). 

[20] X. Ren, I. Leonov, G. Keller, M. Kollar, I. Nekrasov, and 
D. Vollhardt, Phys. Rev. B 74, 195114 (2006). 

[21] F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, 
M. Posternak, A. Yamasaki, and O. Andersen, Phys. 
Rev. B 74, 125120 (2006). 

[22] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 
(1997). 

[23] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 
65, 035109 (2001). 

[24] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Van- 





















23 


derbilt, and N. Marzari, Comp. Phys. Commun. 178, 
685 (2008). 

[25] C.-K. Skylaris, A. A. Mostofi, P. D. Haynes, O. Dieguez, 
and M. C. Payne, Phys. Rev. B 66, 035119 (2002). 

[26] L. He and D. Vanderbilt, Phys. Rev. Lett. 86, 5341 

( 2001 ). 

[27] F. Freimuth, Y. Mokrousov, D. Wortmann, S. Heinze, 
and S. Bliigel, Phys. Rev. B 78, 035120 (2008). 

[28] J. Tun g and G. Guo, Phys. Rev. B 76, 094413 (2007). 

[29] F. Schubert, Y. Mokrousov, P. Ferriani, and S. Heinze, 
Phys. Rev. B 83, 165442 (2011). 

[301 Y. Mokrousov, G. Bihlmayer, and S. Bliigel, Phys. Rev. 
B 72, 045402 (2005). 

[31] See http://www.flapw.de 

[32] Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998). 

[33] B. Hardrat, F. Freimuth, S. Heinze, and Y. Mokrousov, 
Phys. Rev. B 86, 165449 (2012). 

[34] S. Ri and S. Ri, arXiv preprint arXiv: 1407.6824 (2014). 

[35] A. Kokalj, Comput. Mater. Sci. 28, 155 (2003). 

[36] H. Usui, R. Arita, and K. Kuroki, J. Phys.: Condens. 
Matter 21, 064223 (2009). 

[37] M. Shelley and A. A. Mostofi, Europhys. Lett. 94, 67001 

( 2011 ). 

[38] C. Li, A. Freeman, H. Jansen, and C. Fu, Phys. Rev. B 
42, 5433 (1990). 

[39] X. Wang, R. Wu, D.-s. Wang, and A. Freeman, Phys. 


Rev. B 54, 61 (1996). 

[40] L. Bellaiche and D. Vanderbilt, Phys. Rev. B 61, 7877 

( 2000 ). 

[41] K. Seemann, F. Freimuth, H. Zhang, S. Bliigel, 
Y. Mokrousov, D. Biirgler, and C. Schneider, Phys. Rev. 
Lett. 107, 086603 (2011). 

[421 R. Bianco, R. Resta, and I. Souza, Phys. Rev. B 90, 
125153 (2014). 

[43] K. M. Rabe, C. H. Ahn, and J.-M. Triscone, Physics of 
ferroelectrics: a modern perspective , Vol. 105 (Springer 
Berlin/Heidelberg, 2007). 

[44] L. Berger, Phys. Rev. B 54, 9353 (1996). 

[45] J. C. Slonczewski, J. Magn. Magn. Mater. 159, LI 
(1996). 

[46] J. Sun, J. Magn. Magn. Mater. 202, 157 (1999). 

[47] A. Chernyshov, M. Overby, X. Liu, J. K. Furdyna, 
Y. Lyanda-Geller, and L. P. Rokhinson, Nat. Phys. 5, 
656 (2009). 

[48] I. M. Miron, G. Gaudin, S. Auffret, B. Rodmacq, 
A. Schuhl, S. Pizzini, J. Vogel, and P. Gambardella, 
Nat. Mater. 9, 230 (2010). 

[49] I. M. Miron, T. Moore, H. Szambolics, L. D. Buda- 
Prejbeanu, S. Auffret, B. Rodmacq, S. Pizzini, J. Vogel, 
M. Bonfim, A. Schuhl, et al., Nat. Mater. 10, 419 (2011). 

[50] H. Krakauer, M. Posternak, and A. Freeman, Phys. Rev. 
B 19, 1706 (1979). 



